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Hoofdstuk V: Numerieke Lineaire Algebra 


Inleiding 


Wij zullen onze aandacht concentreren rond de volgende proe 

blemen : | ’ 

= Het oplossen van een stelsel van n lineaire vergelijkingen 
met n onbekenden. 

= Het geven van een approximatieve oplossing voor een stel- 
sel van n lineaire vergelijkingen met m onbekenden, m < n 
(kleinste kwadraten problemen). 

= Het bepalen van de eigenwaarden van een matrix. (Soms 
is men alleen geïnteresseerd in de absoluut grootste 
of absoluut kleinste eigenwaarde). 


Regelmatig zullen wij ingaan op stabiliteitsaspecten bij 
de bovenstaande problemen en b.v. nagaan hoe stabiel 
eigenvaarden van matrices en oplossingen van stelsels li- 
neaire vergelijkingen zijn t.o.v. afrondfouten. 


Andere problemen die wij terloóps zullen ontmoeten: 
inversie van natrices, berekening van determinanten, bee 
paling eigenvectoren. 


Lineaire ruimten en normen 

Zij X een n-dimensionale complexe (of reële) lineaire 
ruimte. De elementen van X zullen wij vaak noteren als 
kolomvectoren t.a.v, zekere basis Ejsvees Es d.weze 


he 
k 5 
X : td staat voor de vector x zië Xe; 
EN 
Definitie “ Een norm op X is een afbeelding Il l:X > R, 


met de eigens-happen: 

Ci) Ixil = OS x= 0 

Cii) HA xl = LAI Wx (voor alle A e Lof IR) 
(iii) IIx + yll 7 Ixil + IFyll (voor alle x,y e X) 





(17.3) 


(17.4) 


(17.5) 


(17.6) 


(17,7) 


(17.8) 


(17.9) 


(17,10) 
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B1) Ledere norm op X hoort een limietbegrip als volgt: 
hade keen Y*y d dan als voor alle e>0 een N bestaat zoe 
me 


danig dat ly - Yyxll Se voor alle k > N. 
Zoals bekend bectaan voor e!k tweetal normen || K en |l l' op 
X positieve constanten Yi C2 Yz zodanig dat: 
REX Y1 Wall S Ixil < Ya Ixlt 
Lijgevolg: als lin V‚SY in zekere norr dan geldt dit voor elke 
1:09 
andere norr on “< ook. Alle normen on X geven dus aanleie 
ding tot een zolfca limietbegrin. 


: 
n /o 
Bekende nornen on X ziin de HÖldernormen : Wxll= E EA 
<y 7 i=l 
met x= ( zn Al 
X 
n 
Voor pn = 1: Isxlt,= AA | x,l lineaire norn 
1= temen ER EE NE 
n= 2: shae (Ô alfn' kwadratische of euklidische 


norm 
Via limietovergang ziet men dat voor n= ee: 


xt = max Ix‚l sun norm 
hd Er i men 
1Sif:n 


Opgave. Toon aan 
Wxlia SS Mall, S Mm Ixil, 


Walen S Hella S 4, Ill 
Ga na dat Ge gelijktekens inderdaad kunnen ontreden, 


Definitie. Fen innroduct on X is een afbeeldina (ey) 8 XHKSE 
OR) met de eicenschannen 

(i) x‚x) = 0 » x= 0 

(ii) (x,y) = (‚ZJ 

(iii) (x,‚x) > 6 

(1v) x,y) = A. (x4v) (dus (x,Av) =È (x,y) 

(v) (xtv‚z) = (xX,2) + (y‚z) 


Een eindig=dinensionale lineaire ruimte mct inprodukt heet 


euklidisch in het reële, unitair in het complexe geval. 


NA=-61 


n = 
(17.11) Bij een gegeven basis Cyjreveer EC Zullen wij (x,y)= Z Xs eY4 
= 


i=l 
het standaard inprodukt noemen. Bij het standaard innro= 
dukt is ie basis vanzelf orthonornaal std goe = ê, 
(17.12) Elk inprodukt induceert een norr Ixil = (x‚x) * on cen uni= 
taire ruimte X, 
(17,13) Ga na dat de kwadratische norm geïnduceerd wordt door het 
standaard innrodukt. 


j° 


5 18 Lineaire oneratoren en functionalen. 
mmm mre 


(18.1) Definitie Eon lineaire afbeelding van X in zichzelf heet een 
lineair onerator (LO). 


(18.2) Definitie ton lineair: afhoelding van X * € hect een line= 
air functionaal (LF). 


(18,3) Opgave, fia na dat continuiteit van een afbeelding 
onafhankelij! is van Ae norm (val. (17,3)). 


(16,4) Een norn |l |l or “ induceert normen voor een LO A resp, 
LF L, de zgn, seassoeciecerde normen ‘ie wij ook als || |l 


zullen schrijven: 


ied 
(18.5) Hell = sun Axl = zur HAxll 
x0 Wxll=1 
resp. | Lxl 
(13.6) WL = sun TXT = sun | Lxl 
x0 Wx 1 


(13.7) Stelling. De acassocieerde norm hostaat voor elke LO en 
elk LF. Bovendien nag men in (13.5) en (18.6) "sun" vere 
vangen door “nax", 

Bevijs. we zijn klaar indien tre tonen dat een willekeu= 
rigs LO A on n=din. X continu is. 32 cenheidsbol (Ul xii=l) 
is immers commact, an cen continue functie neemt op een 
compact-um cen maximum (voor IlAxll) aan. (Voor Lr's gaat 
het analoon). ‘‘eaans (16.3) is het voldoende de continu= 
iteit to howijzeon t.a.v. zokere norri on X. “ij kiezen 
hiervoor ‘de sus norm (17,5), nadat we ev basis mn X 
gekozen Auen , 
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Stel dat A op de zojuist gekozen hasis de natrix (a, 3) 
hoeft. 
Dan: 


laxll, = max B: Aij” | < max } lass! max |x‚l S 
zin jel 3 1iSn 


Se. Ixil, 
Hieruit volgt direct de continuïteit. 

(18.73) Dus elke LO en LF is continu. 

(18,8) Opgave. De LO's on X vornen een n2-dimensionalc lineaire 
ruimte, de LP's on X vormen een n-dimensionale lincaire 
ruimte, Ga dit na. Toon aan dat de geassocieerde normen 
normen zijn on de respecticvelijke lineaire ruimten, 


(18.9) Opmerkirng.iliet alle normen on de lineaire ruimte der LO's 
or X zijn geassocieerd aan een vectornorn. 
Voorbeeld: stel de LO A heeft rnatrix (a;j)- 
| EE | 
Zij Hall, = (la, zl 9 ‚ de zan. Frobeniusnorm. 
Dan Will, = /m,‚terwijl voor elko geassocieerde norn | Il =l, 
In plaats van IA, schrijft men ook wel TA . 
(18,10) Definitie, Een norr Ill on de lincaire ruimte der LO's op 
X heet nultinlicatief indien voor alle LO's A,B: 
IAB! SAN. IBU. 


(18.11) Opgave. Geassocieerde normen zijn multiplicatief. 
(Gebruik llAxll S HA. Ixii). 


(18.12) Opgave. De Frobeniusnorm (zie (18.9)) is multiplicatief, 


(18.13) Definitie. Dc abs. waarde van de absoluut grootste eigenwaarde 
van een (complexe) LO A heet de spectraalstraal p(A) van A. 


(18.14) Opgave. Voor een geassocieerde norm |l Il geldt: 
HAI >= PCA) (zie ook (20.6)). 
(18.15) Voor iedere willekeurige nonsinguliere LO T 


geldt PLA) = P(TAT"Ì) 


(19.1) 


(19,7) 


(19.3) 


(19u) 


(19.5) 


(19.6) 


(19,7) 


CEB) 


(19.9) 
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Matrices. 


Een lineaire afbeelding van n-dimensionale X in m-dimensionale 

Y kan bij vaste basis worden voorgesteld door een m x n - matrix. 
In het bijzonder kan een LO vp X worden gerepresenteerd door 

een vierkante n xn - matrix, een LF door een 1 x n = matrix. 

Een vector is op te vatten als n x 1 =- matrix. 

We zullen in het vervolg LO's en LF's identificeren met hun ma- 


trix. 


Zij A een matrix Ca;s)s Ast 5 Taen sie 
D . Lti . . her “ad . KJ 5 
efinitie. De geconjugeerde A van A is de matrix (bij) met 


va Tan 


Definitie. De gespiegelde of getransponeerde AT van A is de 


matrix Ce, sj) met ei3 = âsi: 
Definitie. De geadjungeerde of complex-geconjugeerd gespiegel- 
« . . TT 
de A* van A is de matrix (d,) met d;3 ass: 
Opmerking. Analoge definities gelden voor LF's en vectoren 
(opgevat als n x 1 - matrices). De gespiegelde van een LF is 
een vectop en omgekeerd. Voor het natuurlijk inprodukt 
(wi e= E XiVi schrijft men wel (x,y) =z y*.x, 
1-1 


Á. B; 


Opgave. Ga na dat AÌ - A* , A** = AIT - £ A; AB 
(AB)E = BÎ AT ; (AB)* = B*A*, (ATÌj* = (A*)TÌ 


Opgave. Bij het natuurlijk inprodukt geldt (Ax‚y) = (x, A*y)._ 


Definitie. Len matrix A heet hermiets indien A = A*, unitair 


indien A\Ì = A* 


(Voor een reële ruimte X heten zij symmetrisch resp. orthogonaal) 


Hermietse en unitaire matrices zijn voorbeelden van zgn. normale 
p hatten esta saheseaasd 
matrices: 


Definitie. Een matrix A heet normaal indien A*A = AA*, 


(19,10) 


(13.11) 


€19,12) 


(19.13) 


(19.14) 


(19.15) 


(19,16) 


(19.17) 


NA=6L 
Stelling. Voor elke normale matrix A bestaat een unitaire U zo- 
dat U*A U eendiagonaalmatrix is. 
(Voor een bewijs zie bv. (20.14). 


Opgave. Voor elke A is A*.A hermiets, 


Opgave. Stel er is een unitaire U zodat U*A U een diagonaal- 


matrix is. Toon aan dat A normaal is. 


Opgave. Ga na dat in het geval dat A hermiets is de diagonaal- 
matrix waartoe A unitair transformeerbaar is, reëel is. Geef 
het analogon van (19.12) voor hermietse A, 


Voor matrices in het algemeen bestaat de zgn. Jordan normaalvorm. 
aar vorm 


Definitie. Een Jordan matrix is een vierkante matrix met op de 
diagonaal uitsluitend selijke getallen, op de onmiddellijk rechts 
ervan gelegen nevendiagonaal (als het tenminste geen 
1 Xx 1 - matrix is) enen, en verder nullen „A 1 0 0 
0 A10 
\ 0 0aA1 
\D 00 A, 


Stelling. Bij elke n x n - matrix hoort een niet - singuliere 


"Ear te schrijven is als 


matrix T, zodat T 
aan elkaar geregen Jordanmatrices waar- 

van de diagonalen langs de diagonaal b TE 
van A vallen, terwijl de overige 


matrixelementen nul zijn. 


Opmerkingen. We zullen deze stelling 





niet bewijzen. Zie voor een bewijs 
bv. Perlis p. 161, e.v. De Jordanvorm van A is, afgezien van 

de volgorde van de Jordanmatrices, eenduidig bepaald. De 
karakteristieke veeltermen van de Jordanmatrices zijn delers 
van de karakteristieke veelterm van A; men noemt ze de 
elementairdelers van A. Als alle elementairdelers van A lineair 
zijn en A dus op diaprcnaalpedaante sebracht kan worden, heeft A 


n onafhankelijke eigenvectoren, anders niet. 


(19.19) 


(19.20) 


(19.21) 


(19.22) 
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le en evenzovele verticale scheidingslijnen heeft onderverdeeld 
in vakjes, waarbij de k° horizontale en k° verticale lijn el- 
kaar snijden op de hocf diagonaal: 





Opgave. Laten A en B selijkgepartitioneerde matrices zijn. 
Duidt met A3 en B,; aan het blok op de ie plaats van boven en 
de j° plaats van links in A resp. B. Zij het produkt 

Cz AB weer op dezelfde wijze gepartitioneerd. 


Dan is C33 z L Ask Bij: 
Stelling. lim Ak = 0 « p(A) <1 
k:>ro 
Bewijs. 
>. Voor een geassocieerde norm geldt nAXn > otak) z pak, 


=. Schrijf A als T°Î BT, B op Jordanvorm. Wegens Ak = 171 Bk r 


k 0. Ga na dat het voldoende is te tonen dat 
voor eik Jordankastje C van B ck > 0. 


is dan te tonen B 


Schrijf C = AI + b, D een mutrik met op en onder de hoofd- 

diagonaal nullen, zedat Dt = 0 (l afmeting van C). 

Voor k > 1l seláât dus: 
/ 


k 
ck = akr rijn 


dit gaat naar 0 voor k > co, if 


ON 
ik 


: /k 5 à 
Belaln ent 4 voot (Le he 


en 


, 


Opmerking. Als HAI < 1 in enige multiplicatieve norm, dan geìdt 


trivialerwijze lim ak z= 0. Dat voor lAl| > 1 toch lim ak = 0 
k->eo 


. KE ® o 4 
kan zijn blijkt uit (5 ij 


Opgave. IIA > p(A) voor elke multiplicatieve norm. 
N.B. Dit is een uitbreiding van (18,14), 
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Ale] 


(19.23) Stelling. Ë ar is convergent dan en slechts dan als 


(19.24) 


19,25) 


5 20 


(20,4) 


CZ 


(20.3) 


k=0 - 


k 1. 


lim A* = o. De som is dan (I - A) 
Koo 


Bewijs. Noodzakelijk is triviaal. Voldoende 

Alle eigenwaarden van A hebben modulus < 1 en dus is 

det (I - A} #0. > (I - A) | bestaat. 

Wegens (I + A + A? + …. + ASy(T- A) = I - A*Ì geldt | 
p k +1 nd | „1 ak*1cr hs Ay Ì 

(I + A +A' a. AS) =z (1 - A )(I =- A) z.(I = A) = 

1 


k+1 


en het rechterlid convergeert naar (I = A) “ voor k * oo ll, 


Gevolg. I - A is niet-sinsulier als IIAIl ‘< 1 voor een multipli- 


catieve norm. 


vr 


Opgave. Zij 8 Ak ‚k een maechtreeks met convergentiestraal p. 
k=0 zo 
Als P(A) S p dan convergeert E a Ae Toon dit aan. 
k=C 


Matrixnormen 


Reeds in & 19 werden normen voor LO's ingevoefid. Dit zijn dus 
in feite normen voor matrices (vgl. (18.6)).Diverse van deze 
normen (met name de geassocieerde normen) kunnen op meer expli= 
ciete wijze pekarakter! =.erd worden, hetgeen goed van pas zal 
blijken te komen. We zullen enkele identiteiten afleiden voor 
geassocieerde normen van llöldernormet (17.4) e.v. Js 


Lemma. \lAli, = max Zla,;| = maximale kolom atenluut som. 
1 & 
Bewijs: Ga na: HAx Ui “< max Ela;sl- Ixil, . Als het max. wordt 
j Ì 


aangenomen voor j = k, dan geldt het gelijkteken voor x = ej» 
de k° pasisvector. 


Lemma. |A = max Z|a..l = maximale rij absoluut som. 
1 3 &.J 


Med . | 
Bewijs. Ga na HAxll, < ES gl: lxil… Als het max wordt 
aangenomen voor i = k, dan geldt het pelijkteken voor 


Xx z= (sgn(a Jr oneens gnla,, DT als A reëel is, Hoe is het bij 


kl 
Complexe A? 


(20.4) 


(20.5) 


(20.6) 


(20.7) 


(20,8) 


(20.9) 


NA-67 


Lemma. WlAll, = OREL Hr 


XO, ys0 
Bewijs. 
Ke: vu ) is het natuurlijk inprodukt dat dus de kwadratische 
norm induceert (zie (17 WP. redbruik ongelijkheid van Schwarz, 
en bemerk dat het gelijkteken wordt aangenomen voor y z Ax. 


Ga na: |A, z= IA*,,ll U; A U2 Ile = Alli voor alle unitaire 
U, 22 


Opgave. Ga na dat voor elke normale matrix A (dus i.h.b. voor 
hermietse en unitaire matrices) geldt: Ali, = (A). 

(Gebruik de diagonaliseerbaarheid, 4940) ). Ga wwerzena det 

Og > omge ÀL) Il > AcCA) de oogemwarnanden van, Â. 

Lemma. lll, = o(A*A)À. 

WAxlia = (Ax, Ax)” = (A*A af MEL A*A is hermiets dus heeft op 
geschikte basis diagonaalsvcdtnte. Teks anders: Zg U“ A*AU-D (diagonaal ) 
dan man, (Ax‚Ande mer & Ux, Alx) =— max, (ADx, Aux )= koe COTAPA U xx). max Dxzo 


HAI, HAI HAN, HAU : 1 / 
Lemma, AT, * rar” 5 TATr ï TATE liggen alle tussen pe en Jn. 
1 All, 
ee “= de 
Jn TAN 


De grenzen kunnen aangenomen worden. 


Bewijs. Als opgave. Zie ook, (20,15). 


HAI, 


(Aanwijzing voor A AIG = spoor (A*A) = som der eigenwaar- 
F 


den A*A), 
Opgave. Zij A ='0,3 o,s 
Laat zien dat WAloes HAES, , HAN alle > 1 zijn. 


Toon convergentie van nk naar O aan met (19.21) door een dia- 
gonaalmatrix D te vinden zodat not A DI < 1 in een door u 
zelf te kiezen norm. 
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(20.10) Opgave. Toon aan dat voor willekeurige multiplicatieve norm 


en -1 
WAN A+ Dis AM mits 1 BAT <1. 
1 + || BA “|| 1 = || BA “| 


(20.11) Opgave. Zij A een n x n Jordanmatrix (n > 1) met diagonaalcoëff.). 





1 1 -1 1 
Toon aan: TTATeí S |A Is, Ss TR (| Al > 1) 
SIA ee al > ma 
It + van) Ll - Jatn:) 


Toets deze schattingen door de inverse van A expliciet te be- 


palen en de normen te berekeren. 


(ror) Opgave Toom aam dat Voor Lu wi Melliplicahin Merv 
W0+B* pr < Te aba WBAVet. 
Den AT 


Hunt (A+B)_ A (A+B)BA 
(tors) Hat wrr (2040)3 (A+B) Yo AT (A+B TBA, duo MAB) IE v.C Zo 


C2o.i4) Bourge veer (IJ, Lo) 
UU vra. ol kelner Nn Le. maa 7 olemermohcard 
be a hagst ma. Te UV, U verdun , V brven dueheeho 
ds mak T= UV volgeun a) oek V'U'RUVeW du U'RUEZ, 
Z Gorudrahahe Cbelaugrde Hlellig Vl Phan t Ut matin & Aver 
teru haurs U op de dr LE Orne) 
« . “ 
c)A gel ZZ*-Z°Z, H(A) ed. v.d. parduct mids be vree Z2 
vz + bel + + Vzinl®, vry ZZ à Art |z| * =D Ze 0O, Etc 4 


" ko AA 
(Zo,15) Hua von (20.8) Bepaal du 1,2,09 Lef Arn War ( & Ln Vlev 


4 
(ie) 
Got) Belongs tele, (pokauu obbidang ) Veer ens orla man mai 
beotaas zuilkina waDiae UtoV zodat UTAV=S, Z mùt nagalevt waliix „ 
Beus Nauen V zodat V*A*AV-IE ù Oagemand es da el 
V zijn onderdag olirgenaal , & cia ida Kebe ern (Des « Ger 
En U zodat AV-U3 Ünrass UE belikpat & vereenig vulden t leeden mari) 
(Roiz) Men metw or kelen va V bn ol Ù* rug val Tap. 2echto- tn bike Gilain 
Veche van A bj dn Peieyebina waere og + Merde op VAI = maalDee , 
HAT 4/ mel) s nwd Atataat 


(21.1) 


{21.2) 


(21.3) 


(21.4) 


(21.5) 
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Oplosmethoden voor stelsels van n lineaire vergelijkingen met 


n onbekenden. 


Gausscliminatie. 
zeaussciiminatie 


Het stelsel heeft de vorm 
Ar1Xr + eee. + AinXn = Di n 


Ant X1 + …eeeee PH annXn = bn 


sn luidt dus in matrix notatie’ 


A x= b 
met A = Ca;;) ‚ de coëfficientenmatrix 
ps 
{ bi À 
b di: ‚ een gegeven kolomveetor 
ken , 
Xi 
X -f; | ‚ de gevraagde kolomvector. 
Ken! 


De eliminatie methode van Gauss komt in principe hierop neer: 
verdrijf x1 uit de 2° t/m n° vergelijking door elk van deze 
vergelijkingen te verminderen met een geschikt veelvoud van de 
eerste vergelijking; verdrijf op analoge wijze x2 uit de 

3° t/m n° versclijking mbv. de 2° vergelijking etc. Tenslotte 
heeft men het stelsel teruggebracht tot een stelsel met drie- 
hoeksmatrix, waaruit Xn triviaal oplosbaar is, substitutie 
hiervan in de (n-1)-ste vergelijking levert onmiddellijk 


XxX etc. 


n=1 
Vanzelfsprekend loopt Gauss-eliminatie mis als au = O0. “lgemener, 
als we bij Gauss-eliminatie met AC) aangeven de matrix van het 
stelsel dat ontstaat na eliminatie van Xe (dus Ak A) dan 
loopt het proces mis als afk). 0. Wanneer echter het oorspron= 
kelijk stelsel onafhankelijk is kan niet gelden dat 

an t/m al) alle 0 zijn. Men bedenke hierbij dat 


ee Oy 
a van. de vorm /g 


EED pn 
ie 

Beel Ki 
Nl '. 


Derhalve: Als A niet singulier is kan men tijdens het elimi- 


is en dat det A = det ado), 


é--e 


natie proces de vergelijkingen steeds zo omnummeren dat se * 0. 
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(21.6) Men komt zo tot de volgende algoritme: kijk, alvorens Xie te 
elimineren of B #0 is, en zo niet, zoek een 
k p p R : 
hi ‚ 1 = kt1l,...‚n , die niet 0 is, en verwissel de i-de en 


k-de vergelijking. 


(21.7) De gebruikte elementen an noemt men de pivots (=spillen) en 


het zoeken van een a aie niet 0 is plus het vervolgens 


verwisselen neemt men pivoting. 


(21,8) Opgave, Ga na: AC’ = Iek met 
2) _ CL eit 
a, = a35 als Ì = 1 
U) _ (1) 
a a; a : 
kJ > 
a;s 11 1j als 1 1 
A11 
(21.9) Opgave. Laat AS) = (2 ge * 0. 
Ga na: pri), ed ‚ met 
(ke), (k) n 
az aj ie & pes 
k 
n afs’- dik * Aj als i >k. 
akk 
(21,10) pgave. laat AK = (a Kr ln * 0. 
Ga na: alt), L a de) 
met & \ 
Lis . n 
0 E ei) k & 
\ | 
8, | 
Mk 1 
waarin 2 
WK FT Zik Ci>k) 
Ak 


(21.11)Opgave. Stel dat tijdens het eliminatieproces steeds ac 0. 


U weet ant VEbentiuser, ria vernumeven inderdaad zo 1e mits 


aid 
Ga na: det A= a1ì 


(2) (n) 
2 


Zeeeeeceseseee nn° 


. à 
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$ 22. Practische uitvoering van Gausseliminatie. De hoeveelheid 


rekenwerk. 


(22,1)Bij Gauss=eljininatie wordt het stelsel Aczb stapsgewijs 


(22.2) 


(22.3) 


getransformeerd naar stelsels AT) gap Oer, …, n ‚ waarbij 
men voor kzn op een drichoeksstelsel uitkomt. 

In de praktijk voert men gewoonlijk eerst alle transformaties 
van de matrix A uit en berekent pas later de corresponderende 


transfoismaties van het rochterlid b. (we zullen straks zien 
waarom) . 


Stel dat de elementen van matrix A in de rekenmachine zijn 
opgeslagen in een vierkant array C. 

In feite voert men dan alle operaties uit op dit array, 
waarvan we de elementen aanduiden zullen met: c;‚…. 


1) 
In het begin geldt dus: 
(1) 
Sijs ij * 
na het elimineren van Kg: 
sed 
eek 


etc. 
Ne houden ons voorlopig niet bezig met pivoten, dus nemen 


(k) 
aan dat steeds Ak #0, 


Nu is het niet erg inteëssant in de eerste kolom van C, 
voorzover onder de diagonaal gelegen, al die nullen op te 
schrijven die door het elimineren van xj ontstaan. 

We weten toch wel dat op die plaatsen van „2 nullen staan, 


en kunnen die plaatsen van C wel voor iets anders gebruiken. 


Men plaatst nu in e ‚4 (122, on) het getal 
(1) 
m:,= il 
LL ST 


dri 


i.e. de factor waarmee men de eerste rij moest vermenig- 
vuldigen alvorens deze van de i° rij af te trekken. 
Evenzo plaatst meu in de benedendiagonaalelementen van de 


tweede kolom van C de getallen 
(2) 
Beg 2 As ; k 
12 i2 105 sd.) , 


(22,4) 
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dus de factoren waarmee men de 2° rij van a (2) moet 
vermenigvuldigen om xj uit de overige vergelijkingen 
te elimineren. 

Etc, 


Na afloop ziet C eruit als in de nevenstaande figuur, 
waarbij de u's juist de uu uu 


u u 
elementen van de uiteindelijke muuuu u 
driehoeksmatrix uza{n) zijn, mmuuuu 
en de m's de elementen mis | mmamu u u 
voorstellen, {et de getallen mmmmu u 
Miz kunnen we nu het rechterlid mmmmm u 
van het stelsel gaan transformeren ( vgl. (22.1)) 
(1) 
b z= b 
NT ern en 
d SA ge sdk 
EI nt (23 p 
b; = bi mb, i 23 


* 


(22,5) Het is nu wel duidelijk waarom we eerst alle transformaties op 


(22,6) 


(22,7) 


A en daarna pas die op b hebben uitgevoerd( vBl. (22,1)): 


als er meerdere stelsels zijn met dezelfde A, maar met 
verschillende b, dan hoeven we de operaties op A slechts 
éénmaal uit te voeren. 


In het algemeen zal natuurlijk wel gepivot moeten worden, 
stel b.‚v. dat voor k=3 voor het eerst he, = OQO, dus a {320 
en zij ate 0. 

Dan willen we in A3) 


3 . e e 
maar de vraag is wat we dan met de m's in de 1° en 2° kolom 


2 . LJ 
de 3° en de 5° rl) verwisselen, 


van C moeten doen, 

Ge na dat als we in C de gehele 3° en 5° rij verwisselen 
(dus éòk de m's op die rijen) we precies dezelfde matrix C 
krijgen als wanneer we in de oorspronkelijke matrix A de 
3° en 5° rij verwisseld hadden, 

Als later nog eens al) nul prijkt te zijn, maar a * 0 


5 2 sE 
verwissel dan weer de gehele k°“ en r° rij van C, enzovoorts, 


Opgave. Verifieer de bovenstaande beweringen. 


(22,8) 


(22,9) 


(22,10) 


(22,11) 


(22,12) 
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De na Gauss=eliminatie met pivoting resulterende matrix C 
past aldus bij een rij-gepermuteerde matrix A, en we moeten 
het rechterlid dus aan dezelfde permutatie onderwerpene 


Er geldt nog een opmerkelijke eigenschap voor de uiteindelijk 
resulterende elementen van C. 

Zij U de bovendriehoeksmatrix (i.e. een matrix 

met onder de hoofddiagonaal nullen) die op en boven de 
hoofddiagonaal de u's uit figuur (22,y)heeft, dus uzA{D), 
Zij L de benedendriehoeksmatrix met op de diagonaal louter 
enen en onder de diagonaal het resterende gedeelte van C, 


dus de m's uit figuur (22.4). 


Stelling.Indien bij Gauss=eliminatie geen rijverwisselingen 
plaatsgevonden hebben geldt: 

A = L. U. 
Bewijs. In (21.10) bleek : (T-M, JA ‚ met M de 
matrix die bijna geheel uit nullen bestaat en alleen in de 


ACD, Go) 


k° kolom onder de diagonaal op de plaats (i,k) (i > k) het 
getal Mik heeft. 


AD ee n (1) 
Dan: Uz A han (I Mogi CT MH) Â 
Wegens 

EE ik 
(I=) z= II+ Ue (Ga na) 
geldt ê 
eo (1) _ 8 n) 
Az A Ce (I4M) …os (IM) À 


Ga na dat juist Le(I+M,) .…..CI+M__j). 

Het Gauss=eliminatieproces komt dus neer op decompositie 

van de na het eventueel rijverwisselen (pivoting!) verkregen 
matrix / als produkt van een beneden=- en een bovendrieheeks- 


matrix. 


den spreekt van een L-U-decomnpositie, 


Stel A is nonsingulier en heeft een L=U-decompositie, 
Toon aan dat voor 2 vershillende L-Uedrcomposities 

AzlL U en A! z= L'U' een diagonaalmatrix D bestaat 

zodat 'L=L'D en U' z DU, 

Indien van L of U de diagonaal voorgeschreven is, dan is 


de decompositie derhalve uniek. 


€22,13) 


(22,14) 


(22.15) 


(22.16) 


(22,17) 
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We bezien nog de hoeveelheid rekenwerk die nodig is voor 
het oplossen van een lineair stelsel. Als eenheid nemen 
we de accumulatieve vermenigvaldiging (AV), i.e. een 
vermenigvuldiging gevolgd door een optelling. Ook een 
deling beschouwen we voor wat de hoeveelheid rekenwerk 
betreft als 1 AV. 


Opgave.Verifieer de volgende beweringen 
Berekening van een inprodukt in JR, kost n AV 
Matrix x vector kost n? AV 

Matrix Xx matrix kost n° AV 


Het berekenen van elke m (ziel22,3)kost 1 deling= 1 AV., 


dus in totaal in (n-1) WN. 


ke) on A(K*D) zijn gelijk, het 


(k+1) 


De eerste k rijen van Á 
bepalen van de overige elementen van A uit die van 

ad) als de m eenmaal bekend is kost 1 AV per element. 

ek wordt dit dus (n-1)2 AV, voor nl) (n-2)° AV, etc. 


In totaal vergt bepalins, van Am) dus 


Voor A 


n=1 2 n=1 2 1 
TE (n-k)°z 2 k° z E n{n=1)(n=2) AV 
kz1 kz=1 
en de bepaling van het uiteindelijke array C 
In(n=1) + À n(n-1D(n-2) = F n(n°-1) AV. 
De in (22.4) genoemde operaties op het rechterlid kosten 
(n-1) + (n=-2)t....t 1 z In (n=1) AV. 


(n) pn) 


Het oplossen van het driehoeksstelsel Â vergt nog 


eens in (n+1)AV. 


vrt Rat At oi: 
Totaal vergt Gauss-eliminatie: zn + n° -zn AV. 
dee ® an voor n groot, hetgeen verrassend weinig is in 
vergelijking met een 20 eenvoudige bewerking als matrix 


vermenigvuldiging. 


Heeft men de uiteindelijke C eenmaal verkregen dan 
kost oplossen per rechterlid slechts jn (n=1) +}n(n+1)=n? AV. 
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Dit is dezelfde hoeveelieid werk als vereist voor matrix x 

vector; zodat het ook bij vele rechterleden niet loont 

eerst de matrix te inverteren. Immers, ook als men over 

als beschikt kost het berekenen van ant b nos n° AV, 

(22,13)Opzave. Door in Ax= b voor b achtereenvolsens de basis= 
vectoren Bikes en te kiezen, kan men de inverse van A 
berekenen. Ga na dat bij “ebruik van Gauss=eliminatic dit 
kan seschieden in : n° AV. 


(22.19)Opzave,Ga voor verschillende waarden van n (bv. 10 en 100) 
na wat Je setallen betekenen bij een rekensnelheid van 
bv, 50 microsec/AV. 


(22.20)Opsave. be hoeveelheid rekenwerk bij Gauss=elininatie 
steekt heel runstis af bij de hoeveelheid werk vercist 
voor de resel van Cramer wanneer men de daarin voorkomende 
determinanten m.b.v. de permutatieresel berekent ((n+1)!(n=1) 
vermenigvuldigingen allcen). 


5 23. Varianten van Sauss-elininatie: Doolittle, Crout.Choleski. 
nn nne te ee 


(23.1) Er zijn diverse andere methoden die zich slechts van 
Gauss-eliminatie onderscheiden door de volzorde waarin 
de bewerkinsen worden uit:evoerd, hetseen voor computer= 
sebruik zekere orzanisatorische voordelen kan bieden. 


(23.2) Een heel bekende rekenvol orde is afkomstig van Doolittle. 
Stel dat in nevenstaande fizuur 


de aangegeven elementen u en In uuu uu u 
reeds overeenstemmen net dit van muuuuu 
het array C uit(22.4u), zoals ze mmett+ + 
zijn na het elimincren van Xj en X9-. mmt+t + 
Laten de plusjes nor de oorspronkelijke mm + +++ 
elementen van A zijn (afsezien van even= mm + + + + 
tuele rijverwisselinven). 

de u's zijn dus elementen van Ald), de ontbrekende 


63 de . 
elementen van IN ) zijn zonder meer te berekenen: immers 


3.3 


(23.4) 


(23.5) 
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CEN KAL bn an EIJ KLI 
ans s Aij Car e45 en dis S âgs Ci2°25 Cval.(21.9),(22.2)) 
CAN CLI n 
dus ais = EE C1°1j Ci225° 
N à He ' (3) (3) 
Yaar we berekenen eerst alleen > t/m ann « 
s3 63 


Uit Jie elementen za en a, kunnen we immers de 


pivot zoeken en aansluitend eventueel rijen verwisselen 


(val. (22,6 )). ‘le berekenen la: de resterende elementen 


van de 3° rij var INGE -u dat zijn meteen elementen 
(CH) 


van A Tenslotte delen we C43 t/m c door ez3 en 


63 
nebben dan de nieuwe m's (val. (22,3)). 


Dus: uit de in fig. (23,2) aangegeven elementen 
> ie) > 


van C zoals ze zijn na het uuuuuu 
elimineren van xj en Xs kunnen mu uu uu 
we de in nevenstaande figuur aan= mmu u u u 
gegeven elementen van C zoals deze zijn mmm+ ++ 
na het elimineren van X4 t/m X3 mmm ++ 
bepaien. mmmett + 


det zal duidelijk zijn hoe men het Doolittle-proces 

van begin tot eind doorvoert. 

Het is numeriek equivalent met Gauss=eliminatie volgens 

8 22, d.w.z. de afrondfouten in de berekende coëfficienten 
van L en U zijn bij beide methoden gelijk (mits men bij 
beide methoden steeds dezelfde pivots kiest). 


Het belang van Doolittle's methode is drieledig: 

- Men hoeft niet telkens tussenresultaten in Cop te 
bergen. Als men afziet van rijverwisselingen kont op 
elke plaats van C zelfs maeteen het uiteindelijke 


element terecht. Dit is overigens maar een gering 


voordeel by gebruaik Van ten, Comp. Hat was belanggk bij her pore, mat fel uden 
- De coëfficienten ali) berekent men bij Doolittle in 


feite door van asf’ een inprodukt af te trekken. (vgl. (23.2)). 
Datzelfde zebeurt bij Gauss kennelijk ook, omdat daar 
dezelfde bewerkinsen worden uitgevoerd, maar dat inprodukt 
wordt bij Gauss stapsgewijs opgebouwd (per stap Gauss 


1 term erbij). 


(23.6) 


(23,7) 


(23,8) 


(23,9) 
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Nu is bij de meeste programmeertalen de berekening van 
een expliciet uitgeprogrammneerd inprodukt ( en dat heeft 
men bij Gauss), nogal tijdrovend, 

Vandaar dat men vaak speciale snelle routines heeft voor 
het berekenen van inprodukten, en die kan men bij 
doolittle gebruiken. Dit spaart gewoonlijk een factor 

op de rekentijd. 

Als men bij Gauss in dubbele precisie wenst te rekenen, 
en de uiteindelijke matrix Am) weer op enkele precisie 
wil afronden,heeit mein voos ue hele matrix dubbele precisie 
geheugenplaatsen nodig is, i.e, 2x zoveel ruimteals bij 
enkele precisie. Bij Doolittle kan men het accumuleren 
van de inprodukten in dubbele precisie uitvoeren en de 


(n) 


verkregen elementen van A na afronding in enkelr 


precisie in het geheugen opbergen. 


Als men afziet van het pivot zoeken is Doolittle in feite 
een directe manier om de LU-decompositie van A te bepalen 
(zie (22.11). | 

Men komt dan nl. op precies dezelfde inproduktformules als in 
CAE 2 8 a;j* z 5 lik je 

Bij Doolittle had men de eis dat L een diagonaal van louter 
enen had. 

Een soortgelijk proces is dat van Crout. 

Bij Crout wordt een decompositie L'U' bepaald, waarin L' 
een benedendriehoeksmatrix is, U' een bovendriehoeksmatrix 
met louter enen op de hoofddiagonaal. Ga zelf na hoe de 
coëfficienten van L' en U' in één rekengang bepaald kunnen 
worden en hoe het pivotten dient te geschieden. 


Voor syinmetrische positief definiete matrices, i.e. symmetrische 
matrices met de eigenschap dat (Ax‚x) > 0 voor x * 0, 

gebruikt men veelal de methode van Choleski. 

Deze berust op de volgende stelling: 


Stelling. Zij A symmetrisch, positief definiet. 
Dan bestaat een benedeniriehoeksmatrix G zodat 
T 
A = G.G . 


(23,10) 


(23,11) 
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Bewijs. “fet inductie kan men tonen dat A zonder rijverwis- 
selen een L U-decompositie heeft waarbij de hoofddiagonaal 
van U uit louter positieve eiemnten bestaat. 
Er bestaat eert diagonaalmatrix D( die dan noodz. positieve 
diagonaaleleinenten heeft) en een bovendriehoeksmatrix U' 
met louter enen op de hoofddiasonaal zodat U = D U', 
Dan Az LD U', 
Ga na dat L,‚, D en U' eenduidigbepaald zijn (zie (22,12)). 
Symmetrie impliceert : Li = U', zodat 

Az=LDLF ev pepirfer De. (LDPE, 
Opmerking. In bovenstáand bewijs is Dè gedefinieerd door 
Dè. D* = D. Ga na dat Dè een Jiagonaalmatrix is met 
reële coëfficienten. 
Opgave. Ga na dat de decompositie in (23,9 ) uniek is, op 
het teken der kolommen na. (vgl.(22,12)). 


(23,12) Bij Choleski bepaalt men G = (zij) op de volgende wijze: 


bepaal 11 


bepaal 81 úà > 1) 
bepaal #22 

bepaal Big (i > 2) 
etc. 


C3,13) Merk op dat net vanwege de symnetrie van A overbodig 


is het boven- (of beneden-) diagonaal gedeelte van A 

in het computergeheugen op te slaan. De elementen van 

de ge kolom van G kan men schrijven op de plaatsen van de 
elementen Ais á > j) van de hak kolom van A (Ga na). 

Men ziet dat aldus de hoeveelheid werk en geheugenruimte 
nagenoeg gehalveerd wordt (in vergelijking tot b.v. 
Doolittle en Crout,.)- 

Een tweede voordeel is dat men bij deze methode niet 
behoeft te pivotten, 
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8 24 Perturbatie van stelsels lineaire vergelijkingen. 
eee VeErgerljkingen. 


(24.6), (24.7), (24.8) en het bewijs van (24,9) kunnen door nic* 
hoofdvak-studenten wiskunde worden overgeslagen. 


5 246,1. We zullen later zien, dat de afrondfouten tijdens het oplossen 
van het stelsel vergelijkingen Ax = b gemaakt, geïnterpreteerd 
kunnen worden als een geringe verstoring van de matrix A. D.w.z. 
door berekening verkrijgt men i.p.v. de ware oplossing x een 
vector X waarvoor geldt (A + AA) Xx = b. 
Ook gebeurt het vaak dat de elementen van A en b zelf door be- 
rekening of meting verkregen zijn en dus evt. (licht) verstoord. 
We vragen ons daarom af wat het effect van perturbaties 8A en 
êb is op de oplossing van het stelsel Ax = b. 
Zij dus 


(24.1) (A + 8A) XZ = b + 8b. 
met A + 8A non-singulier. 
We zijn geïnteresseerd in Xx - x= Ax. Er geldt: 
Ax = (I + AT sa) WTsn - sAx). Als nu HATÎgAN << 1 in een re- 
levante geassocieerdenorm dan geeft 


(24,2) öx = A"Ì(&b- 8Ax) 
een goede benadering voor Ax. (ga na) Deze êx zullen we voor- 
taan beschouwen, 


(24.3) Opgave. Ga na, dat voor de differentialen dA, dx, db geldt: 
dx = A°ÌÎ (db - dax). 
Dit betekent, dat 8x een cerste orde benadering van Ax is. 


Uit (24.2) volgt: 
(24.4) ext < lATÍN Heal Ixt + HATÎH HeDh, 


We vragen ons af of het gelijkteken in (24,4) kan optreden, m.n. 
hoe grof ongelijkheid (24.4) is. Daartoe introduceren we het be 
grip maximaliserende vector. 


(24.5) Definitie. We noemen xe X voor een gegeven vector norm maximaliservin 
vector van een LO A of cen LF L als voor de seassocieerdenorm 
geldt laxl = HAN xl pesp. |[Lx| = ALK Ext. 


(24.6) 


(24,7) 


(24,8) 


((24,9) 


8 24.2 


(24,10) 


(24,11) 


(24,12) 
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Stelling. (Hahn - Banâch) Ieder reëel LF F gedefiniëerd op een 
lineaire deelruimte Xo van een reële lineaire ruimte X is met 
behoud van norm uit te breiden tot een LF op X. 


Bewijs. Ieder boek over functionaalanalyse. 


Stelling. Voor een gegeven xeX,xfOen gegeven norm op X is er een 
LF L op X waarvoor x maximaliserende vector is. 
Bewijs. Zij V de verzameling van de veelvouden van x. Definieer 
op V een LF L door Ly z= alxl als y= ax. Dan is [Ly | = ly ! 
voor alle Y & V, dus ÏLia4 en elke y & V is maximaliserende 
vector. Na uitbreiding van L volgens Hahn - Banach geldt nog straits 
ILh = 1, en dus blijft x een maximaliserende vector. 

zee | 
Gevolg. Voor x‚,c € X en een gegeven norm op X is er een LO B op 


X waarvoor Xx maximaliserende vector is en waarvoor Bx = Cc. 


Dit ziet men in door B te definiëren als By = (Ly Je, L het in 
(24,7) bedoelde LF. PN 


Stelling. In ongelijkheid (24,4) kan het gelijkteken optreden 
voor gegeven Ï8Al en Ü8bl en voor een gegeven norm op X en alle 

A en b. 

Bewijs: Zij m een maximaliserende vector van A bm = 1. Kies 
8b = M8Dll.m. Op grond van (24.8) is er een 8A van gegeven norm 
Isa zo, dat x maximaliserende vector van SA is en Ax = =ll SAI xl m. 


A posteriori schatting voor de perturbaties van het oplosproces. 


Vooraf enkele definitees 


Definitie. A = (a; 5) en Bz (b; 5) 
3 $ 
We definiëren A > B indien a, ; > b. ; voor alle i en j. 
1, 1, 
Definitie. a,b vectoren. 


We definiëren a > b indien a, > b; voor alle i. 


Notatie. Zij A de matrix (a, 5) 
bd 
met |A| geven we aan de matrix (la; ;- Zij b de vector (bi); 
b] 


dan is [bl de vector (Ip, 1D. 
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(24,13) We beschouwen en gegeven stelsel vergelijkingen Ax=b. Stel nu, 
dat de a,,; met onzekerheden Aa,; en de b, met onzekerheden Ab. 
_ belast zijn. Men zou dan alle stelsel A x= B met |A- Al < AA 
en |B-b| < Ab mogelijke stelsels kunnen noemen, en hun ennn 
singen mogelijke oplossingen. 
De vraag, die we ons nu stellen is: 
als de vertor y gegeven is, hoe is dan te zien of deze een mo- 


gelijke oplossing is. De volgende stelling doet hierover een uit- 
spraak. 


(24,1) Stelling. Er is een A met [SA] < AA en er is een êb met 
|8b| < Ab terwijl (A + 8A)y = b + êb d.e.s.d. als 
Ir) < AA !y | + Ab met r = Ay = Db. 


bewijs : 
=: (A + 6A)y - (b + ôb) zr + Ay = Db = 0, 
hee pj -Eôa,, «ys tÔb. iz 1(1)n. 
 Iesl <maa, SÌ 1ysl « Leo; 
of ook: [rf <& |&A| ‘In + Ln 
: Írl < Aa |yl + Ab. 
= : [rl < AA |-yl + Ab. Voor de i9°® coördinaat r, betekent 
dit: 
n 
- E Aa, 5 1y3l - Ady Cr, CE Aa, | y;l + Ap, 
j=1 is a: 
Men gaat nu eenvoudig na dat er 6e; „3 en êb; te kiezen zijn 
zo dat £ zl 
|8a; 3! < da; ; en |éb,| < Ab,. 


(24,15) Voor de Häldernormen geldt: 


19. Exp = 1 |x| 1 

2%. als 0 <xS< y dan Ext <t yt. 
3°, HAB S< 1 JA} 1 il 
u, als 0 <AS<B dan BAE, < Isl 
Bewijs dit. 


(24.16) Opmerking. Stelling (24.14) staat toe numerick na te gaan of 
vooraf gekozen AA en Ab voldoen. 


Zijn AA en Ab bekend, dan kan men vervolgens BAAl en lApl 


s 24.3 


(24,17) 


(24,18) 


(24.19) 


(24,20) 


(24,21) 


(24,22) 


lsxh 1  H6AN Usb 
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bepalen in een geassocieerde Hölder norm en daarmee een boven- 
grens voor SAÏ en Mébl . 


Conditie getallen. 


Wij gaan uit van (24,2): êx = A Ìsp = 6A.x). 
afschatten geeft: 


bext < MATIG haxtiePl + Heal Ixt >, 


zodat voor de relatieve fouten geldt 
Tb 


Zij lAt. 1ATÌ4 = C(A), dan geeft C(A) aan hoe de relatieve fouc 
in A en b door kan werken in de relatieve fout in x, 


C(A) noemt men het bij de gebruikte norm passende conditiegetal 


van À. 


Opgave. Toon aan C(A) > 1. Ga na dat C(A) = C(aA) voor a een 
scalair # 0. 


Opgave. Zij A unitair. Toon aan dat t.o.v. de euclidische norm 
geldt C(A) = 1. 


Opgave. Zij A niet singulier. Toon aan dat t.o.v. de euclidische 
norm geldt C(A) = HEX A) : 1 eaten 
minv |A AYT ‘aanw: HAI = inf HA x 1, 
Wxll=1 


We gaan de meetkundige betekenis na van het conditiegetal pas- 
send bij de 2-norm (z= euclidische norm). Zij A; steeds de i=de 
kolom van de matrix A, A; opgevat als heee uit de IRÀ Even- 
zo zij A, de gespiegelde van de ide rij van A dis opgevat als veec- 
tor van de RF 

Dan geldt: 


Stelling. Zij a, de euclidische afstand van A; tot het hypervlax 


opgespannen dese de overige À; en Zij À niet singolier 


1 
Dan is hand = En e 
Â. 


(24.23) 


be As SCA) € Im ) 
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Bewijs: We nemen eenvoudigheidshalve Ì = 1. 


De afstandsvector staat loodrecht op Aa sees sÂn: Wegens 


A A; =0 voor j # 1 is de afstandsvector dus een veelvoud van 


Ar. 





De n 
Stel Ai = AA + E As As. 
| en 
J= 
dan geldt: a; = |Al. HÀr ha. 
Uit Bs “A, z= 1 volgt [Al 1À, 15 ss 
Dus |Af = Rn Oh Ems 
ha, 15 HA, La 


Stelling. Voor alle i en een zekcre j geldt: 


\TATa  , \ rats. KK 
Bewijs C(A) = BA. NÀ. ä 
-1 
/ az \ = HAK KA. 
{ K | 2e i 2 
T LL, 


Nu geldt voor alle i MA, ha < IÀl, en is er een j zodat 


Aj > Mp Ml. 


(24.24) Als A, een hoek 8, maakt met het hypervlak door de overige A, 


(24,25) 


(24,26) 


geldt C(A) > SE ‚ Ga na. De hoek 6, is ook de hoek tussen 
de beelden onder dè matrix A van twee orthogonale vectoren, nl. 
de ee basis-vector en een zekere lineaire combinatie van de 
overige basisvectoren. 

Zij nu 8 de minimale hoek tussen de beelden van twee orthogonale 
vectoren, het minimum genomen over alle paren orthogonale vec= 
toren. 

Dan 


0 = 2x arc cotg ( CLA) ) 

Dit volgt uit de Wielandt-ongelijkheid, geformuleerd en bewezen 
in F.L. Bauer, A.S. Householder: Some inequalities involving the 
euclidean condition of a matrix. Num. Math 2. 1960. 

Een meer elementair bewijs is verwerkt in de volgende opgaven. 


Opgave. Zij A symmetrisch en niet -sindulier 
a. Toon aan dat 8 = hoek (Axo, Ax: ) met 
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= C = B en Xt = e, + e. waarbij e ; 

Ag “max © “min $ Fimax min’ ) Emax’ Pmin 

igenvectoren van A ___ bij A Ned 
SEREN a Ë » Smax 1 “max’ “min bij min 

À = maxi À 

xls maxlacal, 

en : ú le Ï= le : Í 
[Aint minfAA)|. max min 


b. bewijs voor dit geval (24,25). 


(24.27) Opgave. Zij A non-singulier. Xo en x, zijn als in (24,26) ge- 
definiëerd bij de matrix A*A, 
a) Toon aan 8 = hoek (Axo, Ax), @ behorend bij A. 
b) Bewijs (24,25), 


(24.28) We geven nog enige gevolgen van (24.23). 
Noemen we een stelsel vectoren bijna-afhankelijk als één van hen 
een kleine hoek maakt met het vlak door de overige dan impli- 
ceert bijna-afhankelijkheid der A; dus een groot conditiegetal. 


(24,29) Len groot conditiegetal impliceert echter niet bijna- afhankelijk- 
heid (vgl.(24.28) ), zoals men inziet iS de hand van een or- 
thogonale matrix waarin één kolom met 10° is vermenigvuldigd. 


(24,30) Een groot conditiegetal impliceert echter wel bijna-afhankelijk- 
heid als de kolommen van A ongeveer gelijke norm hebben. Laat 
nl. de kolommen van A maar prvcies gelijke norm hebben, dan is 
Inl, < bal, mda, voor elke j zodat uit (24.23) volgt dat 

N n 
voor zekere j geldt C(A) < a 


(24,31) Opmerking. Bij andere normen zal de interpretatie van het con- 
ditiegetal natuurlijk anders zijn. Men bedenke echter dat de 
verhouding van verschillende normen op een eindig dimensionale 
ruimte begrensd i8. Met name voor de 1=, 2- en o=-normen en voor 
de gebruikelijke waarden van n (zeg n < 100) zijn deze grenzen 
< 10 of 100 (zie (20.8) ). Waar het hier slechts om kwalita- 
tieve beschouwingen gaat zal overgang op andere normen dus geen 
essentiële wijzigingen veroorzaken, en zullen factoren als de 
genoemde ons weinig belang inboezemen. 


(24.32) Gewoonlijk zullen schattingen voor Héxl, m.b.v. (24,17) uitgaan- 


(24,33) 
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de van 8A en ii shed 5 orsnschetElngen zijn. Zo geldt b.v. 

voor 8A = AA: x! LE ‚ Ongeacht het conditiegetal van A. 
X 

(24,17) is echter wel scherp in de zin van stelling (24. 9/ Het 

conditiegetal is slechts dan een goede maat voor de mogelijke 

versterking van de relatieve fout als alle &A van gegeven norm 

mogelijk zijn. 


Opgave 
100 0 10 dt À 
A =| 0 100 -u0 b=l0o |. 
10-40 48} 10 | 


In twee decimalen nauwkeurig rekenend vindt men als opl. van 
Ax = b de vector x* 

x-Ï = (-0,88 , 3,6, 8,9). 
Geef een bovengrens voor Ux - x*il, 
Ga na, dat men mag aannemen dat A + 8A weer symmetrisch is en 
nullen heeft waar A nullen heeft. 


/ tav, òA, ruit t.a.v. Âl- 


(25.1) 


(25.2) 


(25.3) 


(25.4) 


(25.5) 
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Numerieke uitvoering van het proces van Gauss - eliminatie. 


In de eerste stap Gauss - climinatie wordt bij de 2° t/m n° rij een 
vector in de richting van de eerste rij opgeteld. Als deze vector 


“groot” is t.o.v. de rij in kwestie (b.v. als Er klein t.o.v, 
ie terwijl de overige elementen van de 1° en i° rij dezelfde groot- 


te orde hebben, (vgl.(21.8) ) dan heeft deze rij na bijtelling 
bijna de richting van de eerste rij m.a.w. in deze situatie wordt 
het stelsel bijna - afhankelijk gemaakt. Aangezien men kan ver- 
wachten dat bijna - afhankelijke stelsels gevoeliger zijn voor 
perturbaties (en daarmee wellicht ook voor afrondfouten) dan niet 
bijna — afhankelijke stelsels is dit ongewenst. (vgl. (24.28) ). 


Voorbeeld van bijtelling "grote vector": 
0,00067 xi + x2 = 1,999 
0,5 x1 + 2x3 = 3 
Rekent men in 4 decimalen (vgl. 5 u) dan luidt het stelsel na 
1 stap Gauss = eliminatie mct a als pivot: 
0,00067 x1 + xa = 1,999 
-Ju4,3X3 = -1489 
zodat x1 = =2,985 
Xa 5 2,001 
De exacte uitkomsten (afgerond op 4 decimalen) zijn echter 
Xi = =2,001 
2,000 


X2 


Opgave. Ga na in (25.2) dat als x2 = 2 + e,‚ uit de eerste vgl. 
volgt: xj = =1, 493 = 1493.e. 


Voorbeeld van een stelsel dat zelf al bijna afhankelijk is: 
Passen we als in (25.2) Gauss = eliminatie toe op het stelsel 
1,01 x1 + 3,5 xa = 7,82 

0,34 xXx, + 1,1 x2 = 2,61 


dan vindt men, rekenend in 4 decimalen, x1 = 3,432, x32 = 1,280. 
De exacte antwoorden (afgerond op 4 decimalen): x1 = 3,326, Xa = 1,314. 
(1) 


Het ligt dus voor de hand als eerste pivot niet zo maar een a 4 
te kiezen, die # 0 is (zie (21.6) ), maar er naar te streven dit 
zo te doen dat de t.g.v. het eliminatieproces bij een rij opgetel- 
de vector niet groot is t.o.v. deze rij. 


(25.6) 


(25.7) 


(25.8) 


(25.9) 


(25.10) 


(25.11) 


(25.12) 
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Dit kan als volgt gerealiseerd worden. 
Vermenigvuldig alle vergclijkingen met zodanige factoren dat in de 


coëfficiëntenmatrix elke rij ongeveer gelijke norm 5) heeft 
(in een willekeurige, maar vast gekozen norm). 
Neem nu als pivot ai met ai) = max lais”! (zie 21,6) en ver- 


k 


wissel de eerste en de i“ vergelijking. 


Ga na dat bij de eliminatie van x, aldus bij elke rij een vector 
wordt opgeteld van hoogstens norm 6. 


Het brengen van de matrixrijen op ongeveer gelijke norm noemt men 
rij - equilibratie. 


Opmerking. Wanneer men in een ongeëquilibreerde matrix het abso- 
luut grootste kolom element als pivot neemt, behoeft dit niet het 
beoogde effect te hebben. (Ga na). 


Bij de eliminatie van x2,... volstaat men met pivoting (dus hier 
het zoeken van de absoluut grootste nd ‚ Ì > k, plus de bijbe- 
horende verwisseling van rijen) en laat men equilibratie achter- 


wege. 


Dat het eliminatie proces ook zonder tussentijdse equilibratie 
"stabiel" is (in die zin dat bij de rijen van de oorspronkelijke 
matrix geen willekeurige grote combinaties van de overige rijen 
worden opgeteld) kan men als volgt inzien. 


Stel dat de vergelijkingen reeds zo gerangschikt staan dat tij= 


dens het proces zonder tussentijdse rijequilibratie 

ed K z ae voor alle k, d.w.z. dat de pivot steeds op 
de juiste plaats verschijnt en geen rijverwisseling nodig is. 
Zij ge de i° rij van Aq). 


ge ontstaat uit ri) door bijtelling van een vector 
Ke Gj) 
B Bee ‚ je.;f <1. (i > k) 


(25.18) 


(25.14) 


(25,18) 


(26.1) 


(26,2) 
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Met inductie ziet men in dat de bijgetelde vector een norm 

S (akr1o4) 8 heeft (vgl. ook (26. 13), hetgeen voor grotere 

waarden van k erg kan oplpen. (Rekenen we b.v. in 13 decimalen 

dan kan blijkbaar voor k > 40 pe geheel ondergaan in wat er bjgehld 
wordt). In elk geval is wit er hij een rij opgeteld wordt begrensd. 


En in werklijkheid zal het natuurlijk zeer zeldzaam zijn dat 


k=1 3 a 
inderdaad U £ age z ck 1_1)e. Maan het kam hel: AS 
j=1 —4-110f 
Dad Ì A41 
… EN EN « =A Arnd 
Overigens is een goede vergelijking tussen het proces met en zon- 
der tussentijdse equilibratie … moeilijk, en nog niet 


uitgevoerd. 


Het hierboven beschreven proces, waarbij men alleen aan het be- 
gin equilibreert en pivot op de aangeduide manier is wel het 
meest gebruikte proces om vergelijkingen op te lossen. Men noent 
het: Gauss - eliminatie met cquilibratie en partial pivoting 
(partial in tegenstelling tot complete, welk geval men niet het 
absoluut grootste element in de aan de beurt zijnde kolom op- 
zoekt maar het absoluut grootste van de hele nog te bewerken 
rechtsondermatrix).. 


Het effekt van afrondfouten 


Het door berekening verkregen resultaat van een proces is vaak 

te interpreteren als het exaete resultaat behorend bij enigzins 
verstoorde beginwaarden (vgl. 524.1). Dit is ook het geval bij 

het Gauss - eliminatieproces met partial pivoting. Een uitvoe- 

rige analyse is gemaakt door Wilkinson die voor dit proces, 


evenwel ongeacht of de matrix geëquilibreerd is, heeft aange- 
toond: 


Stelling. De via Gauss - eliminatie met pertial pivoting bereken- 
de oplossing van het stelsel Ax = b voldoet exact aan: 
(A + AA)y = D 


met « (Aat i E 
PN < 1,01. &(n). glA)E 

















(26.3) 


(26.4) 


(26,5) 


(26.6) 


(26.7) 


(26.8) 


(26.9) 
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Waarin ò(n) < n? + 3n? 
VAN = max la;;! 
max las“ 
gla) = i,j,k,' 2) 
bal 


De grootheid g(A) geeft dus aan hoeveel maal zo groot het groot- 
ste element te eniger tijd in het eliminatieproces is als het 
grootste element van de oorspronkelijke matrix; g(A) is dus een 
soort groeifactor. Merk op dat g(A) > 1 wegens aj = as) 


Uit (24.4) volgt nu: 


ile < tatin. haat 


1 f _ 
5 < lA kb WAL. 1,01. dln). g(A). E 


kasd 
oo 


De grootheid c'(A) z HAT. VAR is weer een soort conditie 
getal. 


Opgave. Ga na C'(A) > 1, C'(AA) = C'(A), en 
1 C'(A) C' (B) 
C'(AB) > = ax n re |} 
a creBD) * era d) 
(Hint: bedenk dat nlAl, cen multiplicatieve matrix norm is.) 


We bezien de invloed van rijequilibratie op C'(A). Ga na dat 
rijequilibratie van een matrix neerkomt op linksvermenigvuldi- 
ging met een geschikte diagonaalmatrix, 


Stelling. C'(DA) is minimaal als in DA alle rijen gelijke 

ee « normen hebben. 

Bewijs. 

Stel dat in A reeds alle rijen gelijke © « norm 8 hebben. Dan 
is IDAI= ID. 6 = Ipt tal 


m 

Wegens HA “Îp Inlog la ÍpÎk pb geldt: 

C'(DA)-= HA"ÍDTÎN . IDAN_ > HA Îea lat. "Pim z= C'(A). 
5 oe m m TDF 


oe 


(26,10) 


(26,11) 


(26.12) 


(26.13) 
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Derhalve: equilibratie minimaliseert de factor C'(A) in het ú 
rechterlid van (26,5). Nemen we aan dat A al geëquilibreerd is 
dan kunnen we zeggen dat geen enkele ‘andere rijschaling een 
rechterlid voor (26.5) kan opleveren dat meer dan een factor 

g(A) kleiner is dan voor A zelf (immers g > 1 voor elke matrix). 


De vraag is nu nog wat eaquilibratie met g(A) doet. We merken op 
dat t.g.v. partial pivcting g(A) < ghel, Om dit in te zien be- 
schouwen we het elimineren van xj. Stel dat |a,,! = max lag | 

k 


(zie (25,5) ) zodat men ayals pivot gebruikt. Dan wordt bij 
elke volgende rij hoogstens de 1° rij opgeteld of afgetrokken, 


zodat max ja; add S 2 max |a,;l. 
ij 
i,j isj 


Analoog bij de volgende eliminatie stappen. Men kan zelfs een 
matrix aangeven waarbij zn wordt gehaald. Men zal echter 
begrijpen dat dit een uiterst ongelukkige samenloop van omstan- 
digheden vercist, en gelukkig blijkt in de praktijk slechts hooest- 


zelden g(A) > bh, of de matrix geëquilibreerd is of niet (het 


"wonder van Wilkinson": „Wilkinson heeft dit ontdekt) 


Men kan derhalve zeggen dat vquilibratie vrijwel steeds het rech- 
terlid van (26.5) tot zijn minimum veduccert op hoogstens een 
factor hk na. Dit is een indicatie voor het praktische succes var 
de methode. 


Merk op dat men tijdens de berekening kan vaststellen wat g(A) 

is. Mocht nen g(A) te groot vinden, dan kan men overstappen op 

andere opiosmethoden. Deze, alsmede vele nadere details van het 
voorafgaande vallen echter buiten het bestek van dit college. 


(27.1) 


(27.2) 


(27,3) 


(27.4) 


(27.5) 


(27.6) 


(27,7) 
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Naverfijning 


B 5 .e 

Met xj e.d. (is1,2,...) zullen we hier het i =element van eer 
. . d e e. . 

rij van vectoren aanduiden en niet een coördinaat. 


Stel dat we via Gauss - eliminatie met partial pivoting voor 
Ax = b een benaderde oplossing xj vinden. 
Enige verbetering van deze waarde is soms mogelijk door nauwkeu=- 
rige berekening van het residu rm, = b = Ax: en als nieuwe bena- 
dering te nemen: 

X3 = Xi + Va 
waarin v, de berekende oplossing is van Az = r: . Bedenk dat 
Z = vi nu snel bepaald kan worden omdat de LU - decomp. van het 
stelsel na eenmaal oplossen van Ax = b bekend is (zie ook (22.45) 
Het proces kan men herhalen met x2 etc. 


Er wordt zo een rij Xj , X2 » X3 seees Xjes Xe? et VAN benaderde 
oplossingen geconstrueerd waarbij Xie 5 *K * Yi OP Vi de met 
Gauss = eliminatie berekende oplossings vector is van het 

stelsel Az z= Pk 5 b = Axy 


‘De vraag is nu: convergeert de rij x, naar de exacte oplossing 


van het oorspronkelijke stelsel? 


Wegens (26.2) kunnen we aannemen dat Vi exact voldoet aan een 


vergelijking: (A + (AAD) Vi 7 Fr 


Stellen we B, = A” (a) aan geldt: 


ACI + BĲ) Vi * Fk 
Stelling. Laat in (27.6) voor alle k B, < e < } in zekere vast- 
gekozen geassocieerde norm. Dan convergeert de rij X, naar de 


exacte oplossing van het stelsel Ax = b. 
1 


Bewijs Zij a exacte wortel van Ax = b, d.w.z. Gz A “b, Laat 
he 2 - X 

a oe Mnlar dke : 

had ACI + BĲ) %Xeag * ACI + BĲ) X + Fy * AB, Xx * D 

- (1 + BĲ) Krat * Bik * & 


® CI B) h B, h 


k k 


k+1 


(27.8) 


(27.9) 
(27.10) 


NA-92 


Wegens LB! <1 is CI + B) nonsingulier({19.24)) zodat: 


se „} 
hia4 = (I + BĲ) . B h 


> In/l < HCT + BOTÎN. HBI. An U 
_k+1 k k k 


Met (20,10) volgt: 


IB, & 


k 


zodat voor ip! SS e < } geldt: 
€ 
Inr44! < TE In! < Int. 


Corollarium. Naverfijning bij Gauss = eliminatie met partial pi- 
voting convergeert als 

1,01. C'(A). dln). glA). E <3 
Bewijs. (26.2), (27.5), (27.6) en (27.7). 


Opmerkingen. 


Er is aangenomen dat arithmetische fouten alleen optreden tijdens 
het oplosproces en PEN X + Vi ‘exact! berekend worden. 

Rekenen in enkelvoudige precisie is dan niet toereikend, De af- 
rondfouten bij de berekening van inprodukten zullen de evaluatie- 
fout in r, =b - Ax, (relatief) aanzienlijk doen oplopen. Men 

mag niet meer hopen dat de oplossing van dit stelsel Az = Pe (met 
Gebruikt men evenwel een dubbellengte - procedure voor berekening 
van inprodukten, dan krijgt men Pr voldoende nauwkeurig (althans 
voorlopig). 

Het berekende residu dient als rechterlid in de vergelijking 

Ax = rj, waarvan de oplossing vj behept zal zijn met cen rel.fout 
>= K. (A) E (kr zekere constante). (vgl. (26.5)). 

Mits aan de convergentie conditie van het proces is voldaan zul- 
len de benaderde oplossingen Xx Aldus aanvankelijk op ven adequate: 
wijze door de Vi worden gekorrigeerd. . 

In het gunstige geval trecdt convergentie op d.w.z. zullen vanaf 
zeker moment de decimalen van Xx niet meer gewijzigd worden. Gr - 
tere nauwkeurigheid is dan alleen mogelijk als ook de addities 
Oe + vi) in dubbele precisie worden uitgevoerd. 

Het kan echter ook gebeuren dat na aanvankelijke verbeterinp van 


en ta Geen pag 93 


(27.11) 
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de benaderde oplossingen verdere correcties een oscillatie dur Xi 
doen ontstaan. Het is niet te verwachten dat optellen in dubbele 
lengte dan nog noemenswaardige verbetering brengt. Een fcuten- 
analyse voor het naverfijningsproces vindt men b.v. in C.B. Moler- 


JACM 14 (1967) 316 - 321 


Uit het bewijs van (27.7) volgt nog een indruk van de converg n- 

tiesnelheid. Per iteratiestap vermindert de fout met een faetor 
€ 

ek Ti Ì 

Ook hier is men kennelijk gebaat bij cquilibratie van de matrix, 

omdat dan C'(A) en daarmee dus in zekere zin de convergentiefac- 


tor bij naverfijning (zie (27.8)) geminimaliseerd wordt. 


(28.1) 


(28.2) 
(28.3) 


(28,4) 


… (28.5) 


(28.6) 


NA-95 
Iteratieve methoden. 


Iteratieve methoden (vgl. S11) ter oplossing van Ax = b berusten 


veelal op een geschikte splitsing A = Ai + A3 waarna men het 
stelsel x = (lx) = Aj *b - A7 “Ax door successieve substitutie tracht 


op te lossen. Het is daarbij van belang dat Aj* zonder al te veel 
moeite bepaald kan worden. 


Voorbeelden. 


dezelfde behedonsenoek als A. 

Opgave. Toon aan dat de rij x1 ‚ys ,... met Kier z Aj *b - Ai Apx, 
convergeert dank de wortel van Ax z b indien in zekere geassoci- 
eerde norm HA, ÎA Al <1. 


Men kan aantonen dat Gauss - Jacobi en Gauss = Seidel convergeren 
voor diagonaaldominante matrices, Gauss - Seidel ook voor symme- 
trische positief - definiete matrices. 

Een bezwaar van de iteratieve methoden is dat zij doorgaans slechts 
langzaam convergeren. 

Later in dit college zullen we nog wel nadere bijzonderheden over 
deze processen ontmoeten. Men zie ook b.v. Varga (ch. 3). 
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(30.1) 


(30.2) 


(30.3) 


(30.4) 
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Een belangrijk probleem in de lineaire algebra is het "zo goed 
mogelijk oplossen" van een overbepaald stelsel Ax zb, d.i. een 
stelsel waarin A een m Xx n matrix is, m > n (een dergelijk 
stelsel zal immers i.h.a. geen echte oplossing hebben) Dit pro= 
bleem ontmoet men bijv. bij het uitwerken van de resultaten van 
experimenteel onderzoek, waar men voor een aantal onbekende 
grootheden Xi »e-X, en groot aantal lineaire relaties 

Jai”; z b,, i s= 1,2,..,m kan bepalen, waarvan de coëfficiënten 
echter niet erg nauwkeurig zijn. Of men wenst een functie, die 
men in m punten kent, te vervangen door een (hoogstens) (n-1)-ste 
graads polynoom. 

Onder de “zo goed mogelijke oplossing" verstaan we nu die x waar-= 
voor ÎAx - bl minimaal is in een of andere norm. Met de 2-norm 
noemt men dit een kleinste kwadraten oplossing, omdat voor deze 

x de uitdrukking Tl Ea,jx xj 7 b, Al geminimaliseerd wordt, met de 

ww - norm een Chebichev” - aplossine: In de volgende paragrafen be- 
perken we ons uitsluitend tot kleinste kwadraten oplossingen. 
Bovendien veronderstellen we dat de lineaire ruimten reëel zijn. 


Normaalvergelijkingen (van Gauss). 


De klassieke aanpak van dit probleem is met de zgn. normaalver- 
gelijkingen van Gauss. Dat gaat als volgt: 


Zij A een mx n matrix, m > n, met rang (A) zn. 


Definitie. x € IRF heet een kleinste kwadraten oplossing van 


Ax = b indien voor alle y € IR: lAx - pl, < lAy - bl,. 


Nu beeldt A de RP" af op een n-dimensionale lineaire deelruimte 
van de RÀ die wordt opgespannen door de kolommen (als vectoren 

beschouwd) van A. Blijkbaar moet Ax dan die vector in AIR" zijn, 

die de kleinste 2-afstand tot 
b heeft, d.i. de projectie b! 
van b op ART 

Dus (AR”, Ax - Db) z 0, 

of (RRP, A*Ax - A°bD) z 0. 

We hebben nu bewezen: 
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(30.5) Stelling: x is een kleinste kwadraten oplossing van Ax z b 
als A*A x z A°*b, 


(30.6) Opgave. Bewijs ook rechtstreeks dat lAy - bl, > lAx - bl, 
als A*Ax z A*b en y # x. 


(30.7) Het lineaire stelsel vergelijkingen A*Ax z= A*b noemt men de 
normaal vergelijkingen van Gauss. 


(30.8) Opgave. Als A rang n heeft is A*A positief definiet en dus niet 
singulier. In dat geval kan men dus spreken van de kleinste kwa- 
draten oplossing. Als A rang < n heeft is de kleinste kwadraten 
oplossing niet meer eenduidig bepaald; echter is Ax voor alle 
kleinste kwadraten oplossingen hetzelfde. Ga dit alles na. 


(30.9) Het lineaire stelsel A*Ax z= A*b kan men in het niet singuliere 
geval wegens (30.8) oplossen m.b.v. Choleskí. 


(30.10) De aanpak met Normaalvergelijkingen werkt bevredigend voor niet 
te grote waarden van n (rg n &< 5), hetgeen ìn de tijd van het re- 
kenen met de hand steeds het geval was. Met de komst van de com= 
puter werden veel grotere waarden van n mogelijk en toen bleken 
zeer onaangename verschijnselen. We laten dit eerst zien voor 

het (veelvoorkomende) geval van polynoomaanpassing, en tonen hoe 
het wél moet. Later zullen we dit ook voor het algemenere pro- 
bleem doen. 


6 31 Polynoomaanpassingen. 


(31.1) Op [0,1] is een functie geta p uleerd gegeven in de punten ts vet, á 
homogeen verspreid op het segment [0,1]. Gevraagd wordt een (n=-1) 
ste graads kleinste kwadraten polynoom aanpassing, n=-1 < m. 


(31.2) Schrijf het gevraagde polynoom als: 


(31.3) a, +a xt... + he waarin ag,‚...,ä 


o 1 4 de onbekenden zijn. 


n= 


(31.4) 


(31.5) 


(31.6) 


(31.7) 


(31.8) 


(31.9) 


5 32 


(32,1) 
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m4 , k 
Met A = je AE Za en Xx z /ao hebben 
ke en a 
n=1 1 ij 
1 Erger: De an-y 


we dus een lineair kleinste kwadraten probleem Ax = b, waarin Db 
de m=- dimensionale vector van functiewaarden is. 


De matrix coëfficiënten van de matrix A*A worden dan gegeven 
door 


Tm ols le 
E (t)1°J72, (Ga na) 


id s 
(A AD 3 Ef 


1 ses 
Voor m groot zal gelden: (ATA); 3 ® m { t°) Zat z Tae: Ga na 


Gevolg: voor grote m zal de matrix A*A sterk gelijken op de matrix 
E 


1 | 1 
1 5 be: hele kank 
| 


1 7 
3 .eeeesee n+1 


mx | 
| 
\ 
\ 


5 jet . e 


1 4 
Cr NI Andi Akke zet / 
Zonder de factor m is dit een eindig segment van de Hilbert-matrix. 


(De Hilbert-matrix is de oneindig voortlopende matrix). 

Men kan aantonen, dat het conditiegetal van het n xn - segment 
oe 2°n is. Binair rekenend met 40 - bits kan alleen al het afron- 
den van de coëfficiënten van A*A vanaf n = 9 de oplossing onher- 
kenbaar verminken. 


De juiste oplossingswijze voor dit soort kleinste kwadraten pro= 
blemen loopt over orthogonale polynomia. We geven hiervan een korte 
schets. Voor details zie TA $8 e.v. 


Orthogonale polynomia. 


Bij het oplossen van het in $31 gestelde probleem hebben we geeist, 
dat de oplossing de gedaante p 
n= 


$ . $ a a & + a Xx 


had. 
We kunnen, algemener, eisen dat de oplossing de gedaante 


CI2 2) 


CAR) 


(32.4) 


(32.5) 


(32.6) 


(32,7) 


(32.8) 


(32.9) 


(32.10) 


(32.11) 


(32.12) 


(32.13) 
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dp Pa 4 waa & An-1Pn-1C) 


heeft met P; een polynoom van precies graad i. Geheel analoog aan 


(31.4) kan men weer de formulering als kleinste kwadraten probleem 
geven. Ga na dat nu NGR Pn-1 Ct, ) \, waarbij Eje. tt 


eee ee ese see es ee s e e 


Pot: Pn-1 tr ) 
de gegeven steunpunten zijn. 


We zullen trachten de polynomen P; zo te bepalen dat A*A een pret- 
tige gedaante heeft. 


De matrix coëfficiënten van A*A worden 


). 


k 


m 
Ld 
(APA) js E Pitt) ps Ct 


Een matrix die beslist niet aan het in 531 gesignaleerde euvel 
mank gaat is vanzelfsprekend de eenheidsmatrix. Dat betekent dus 
dat de polynomia {p;} voldoen aan: 


Nu vormen de polynoom van graad aas een lineaire ruimte Pr 
(Ga na!). 


Opgave. Zij voor i = 0,...‚m=1 P; een polynoom van precies graad 
i. Toon aan dat de (p‚}een basis vormen van Pe" 


Voor m gegeven punten Er see est en q1, qs € P me, definiëert 


m 
Cq: ,q2) = B an Ct) qa lt) 


een inprodukt in de ruimte Pm=1 
Opgave. Toon dit aan. 


Blijkbaar impliceert nu onze wens A*A = I dat 
(ps »P; ) = 1 
(pj>Pj) = 0 BE 0Si,j Smi. 
en Pe pracee de aad i fact 
Indien (32.12) geldt’zegt men dat de polynomen Pj een orthogonaal 
stelsel vormen t.o.v. het gebruikte inprodukt. 
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(32.14) Uit de lineaire algebra is bekend hoe men in een lineaire ruimte 
met inprodukt uit een gegeven basis een orthogonale basis con- 
strueert. (het proces van Gramm - Schmidt). Voor de polynomen Pi 
ligt de situatie eenvoudiger, zoals men ziet in TA 58 e.v. 


(32.15) We kunnen de oplossing van ons probleem (het (n=-1) ste graads 
polynoom vinden dat een in m punten gegeven functie in de zin 
van kleinste kwadraten benadert) nu expliciet aangegeven. 

We beschouwen daartoe de m gegeven functiewaarden in t,,...t 


m 


afkomstig van een polynoom f uit P (waarom kan dit?). 


m1 
Dan geldt voor de kl. kw. oplossing Xx: 


(32.16) xs ( (posf) , (pf), ………, (pr-rf) Ne 
en het bijbehorende kl. kw. polynoom is 


(32.17) (po,f) po + (pi sf)p, + «…… + (Prjvf) Pn-4* 


(32.18) Opgave. Voor n = m is de kl.kw. opl. met het Lagrange - interpc- 
latie polynoom. 


(32.19) Wie dit practisch wil toepassen wordt aangeraden TA 58 e.v. te 
raadplegen, waar ook uiteengezet wordt hoe men deze zaken han- 
dig programmeert (met name het gebruikmaken van een drieterms- 
recursie blijkt zeer profijtelijk). Overigens zijn er standaard- 
procedures, die volgens deze lijnen werken, beschikbaar. 
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Ca3s 1) 


(33.2) 


(33.3) 


(33.4) 


(33.5) 


(33.6) 


(33.7) 


(33.8) 


(33.9) 
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De methode van Householder. 


In deze paragraaf zullen we een methode bekijken ter verkrijging 
van een kleinste kwadraten oplossing van het overbepaalde lineaire 
stelsel Ax = b die 

stabieler is dan de methode van de normaalvergelijkingen. 

Het proces berust op het vinden van een orthogonale matrix Q 
zodat QA = U, met U een bovendriehoeksmatrix. 


Allereerst tonen we aan dat zo'n decompositie bestaat. We geven 
een constructief bewijs waarbij de kolommen van de matrix Q ach- 
tereenvolgens worden berekend. 

Centraal in dit bewijs staat de Householder - transformatie. 


Definitie. Voor elke v € R", v #* 0, definiëren we de lineaire 
afbeelding 


T 

H‚=I-2 2, 
Ivt, 

We noemen deze afbeelding de bij v behorende Householder trans- 


formatie. 
Bewijs nu zelf de eigenschappen (33.5) t/m (33.8) 


Hv = Hy voor elke constante a. 


Een Householder transformatie is orthogonaal en symmetrisch. 
Voor uli v =- Hu = U. 
Voor ul v @ H(utv) zu = v. Dus Hy ùo etn Spègelng aa vlak olser OL v. 


Bij elke a en bS mR, a*0,b#*o0 is er een ve IR zodat Ha 
de richting heeft van b, d.w.z. er is een constante À zodat 

Ha = Ab en wel voldoet v = a = Xb met À z + lal, / Ipla. 

Bewijs: Wegens(33.6) zal moeten gelden Az + Tal, /Iblsomdat, van- 
wege de orthogonaliteit In al, = lal,. Wegens (33.8) is het vol- 
doende orthogonale vectoren u en v aan te wijzen zodat a = u + v, 
Ab = u * v. Dit is equivalent met v = at en uz ee 
Hiervoor geldt inderdaad 


(v‚u) = B (a,a) - AF(D,b) } = 0. 


(33.10) 


(33.11) 
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od, a #A6Cmaan dan veldsch elke viv) 


Dus H‚ met v = (a - Ab)/2 voldoet/ Maar dan ook v = a - Ab, aan- 
gezien H, alleen afhangt van de richting van v (wegens (33.5)). 


Op de volgende stelling is de buvervedveïide úecompositie gebaseerd. 


Stelling. Zij a = (as,aa2,...; a)” de eerste kolom van de lin. 
afbeelding A : R'> RP. Zij v de vector 

v = (a: * lal,, a2,..., ah 
Dan is HA een matrix waarvan de eerste kolom uit nullen bestaat, 
uitgezonderd het diagonaalelement. 
Bewijs: Volgens (33.9) voldoet v = a + lalse,, e, eerste basis- 
vector. 


Gevolg: Er is een orthogonale matrix Q en een mxn bovendriehoeks- 
matrix U zodat QA=U. 


Bewijs: Stel dat de matrix AC) 


reeds van de gedaante 


Ks an: DON [ 
| k 
ye) va | A 


oen mot hale an va el ed 
Ck) _ | 


1 

| 

\ m=k 
© wo) 
een kxk bovendriehoeks matrix. Men kan een House- 


holder transformatie nú) toepassen op wei) 


(k) 


is, met ut 


die de eerste kolom 
van W uitgezonderd het diagonaal element nul maakt wegens 
(33.10). 

De matrix 


Jk | Mzk, 
(kx) I | 


S pe nn an nf a tn 


HC [ex 
/ 


IE 
AR 


voert A0) over in alkti), Met de vi e 
sn gtn-2) s{2s 10) pn U zi) 


holdertransformatie uit (33,19 is. 


Of anders geschreven : (n) 
u 
. GA „(e) 
sd 
waarin Q als produkt van orthogonale mxm matrices een orthogonale 


mxm matrix is. 


aC0) z A krijgen we zo: 


waarin no) de House- 


(44.42) 


(33.13) 


(33.14) 
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Omdat Q orthogonaal is dus afstanden 
en hoeken behoudt, is de kl.kw. opiossing van Ax z b ook kl.kw. 
oplossing van Q Ax = Q b , dus van 


(n) 
Sk Qb ‚,QvEe RmrÀ 


(ga na!) 


Een kl.kw. oplossing is gedefiniëerd als die x € Re waarvoor 


LD) el. 


minimaal is. Dit betekent dat x oplossing is van het lineaire 





stelsel 


(n) je 
U X = Qb 


waarbij Q de matrix is bestaande uit de eerste n rijen van de 

8 ' (m) ID 
matrix Q. (Ga na!). Als tamg(A)en a UP mon stonguban « Abs rang (A)< n dan 
ie op zeker momint Tijdone Bat perors da 1 keker van WÓÒ mal „Laar os lest mat an … 
Opgave. Hoeveel AV zijn i.h.a. nodig om un) 
Hoeveel worteltrekkingen? 


Vergelijk dit met het aantal AV nodig om A*A en A*b te berekenen: „/ 


, 


en , te berekenen? , 
1 
' 








p” 
Zy A naïef (complix nd reden ordotknda. qeurven weg 
We leggen nog enige verbandenN Zij weer is edel ds ordner E 
(n) Cn) 
Ks U st A= q' U latan O. 
0 0 


Er geldt A*A = U*U (ga na), maar ook A*A = C*C, waar C de 

nxn bovendriehoeksmatrix is volgens de Choleski=decompositie. Nu 
is C op linksvermenigvuldiging met een matrix S na, uniek bepaald, 
waarbij S een diagonaalmatrix is met +1 of =1 als diagonaalele- 
ment. D.w.z. elke rij van C mag willekeurig met +1 of =-1 vermenig- 
vuldigd worden. C is op linksvermenigvuldiging met S na uniek be- 
paald door de eisen : C*C z A*A en C bovendriehoeks. 

U voldoet hieraan, dus is U een decompositiematrix volgens Choleski 
van A*A en op linksvermenigvuldiging met een S na uniek. 

Voor de eerste n kolommen van de matrix q' voigt dan, dat zij op 
rechtsvermenigvuldiging met een matrix S na uniek bepaald zijn 
d.w.z. de eerste n kolommen opgevat als vectoren in de R” mogen 
naar willekeur met +1 of -1 vermenigvuldigd worden (zelfde of 
tegengestelde richting). Nu vormen deze kolommen een basis van 

de n-dimensionale lineaire deelruimte, opgespannen door de kolom- 
men van A. Bovendien worden de kolommen van A uit de eerste n 


(33.15) 


(33,16) 


(33.17) 
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kolommen van QT verkregen door rechtsvermenigvuldiging met een nxn 
bovendriehoeksmatrix U. Maar dit is ook net het geval als men de 
eerste n kolommen van QF vervangt door de orthogonale basisvectoren 
die het proces van Gramm-Schmidt oplevert, uitgaande van de kolom- 
men van A. Derhalve zijn de eerste n kolommen van Top een fac- 

tor + 1 na identiek met de basisvectoren volgens Gramm-Schmidt. 

Bij de oplossing van het probleem van de nde. graade polynoom aan- 
passing (8 32) kwamen we uit op een orthonormale basis van polynomen, 
op het teken na uniek bepaald en identiek met de basis die het 
proces van Gramm-Schmidt oplevert, Door de identificatie van een 
polynoom met de vector met ais coördinaten de waarde van dit po- 
lynoom in een m-tal gegeven punten tE kan men de n orthonormale 
polynomen gelijkstellen aan de eerste n kolommen van de matrix Q 


Bij de numerieke uitvoering van een Householder transformatie 
heeft men nog de vrijheid het teken in a, + lal, te kiezen 

(zie (33.10)). Zuiver wiskundig gesproken is de tekenkeuze niet 
van belang. Er wordt echter in eindige precisie gerekend. Dan is 
het om numerieke stabiliteit te garanderen noodzakelijk het teken 
gelijk aan sign(a;) te kiezen. Voor een bewijs zie Wilkinson, 
blz. 154 e.v. 


De methode van de normaalvergelijkingen berekende de kl.kw. oplos- 
sing x door het lineaire stelsel A*Ax = A*b op te lossen, Pertur- 
baties in A en b zullen dus maximaal een factor C(A*A) (het con- 
ditiegetal van A*A, in de eucl.norm) vergroot in x doorwerken. 
Uitgaande van de deoompositie verkrijgt men de kl.kw. oplossing 

als oplossing van Ux = Q.b. Men zou verwachten dat perturbaties 
in A en b nu maximaal met een factor C(U) z JC(A*A) (ga na:) ver- 
groot in x doorwerken. Helaas is dit niet helemaal waar. In het geval 
van een groot residu (d.w.z. lAx-bl vergelijkbaar met Hbl) kan men 
aantonen dat weer de factor C(A*A) optreedt. Bij klein residu blij- 
ken de resultaten verkregen met een decompositie gebruikmakend van 
Householder transformaties echter veel nauwkeuriger te zijn dan die 
van de normaalvergelijkingen. Voor details en literatuurlijst zie 
G.H. Golub: Matrix Decompositions and Statistical Calculations, 
Tech. Rep no CS 124, 1969, Stanford University. 

In (33.14) hebben we in feite aangetoond dat de transformatie 

QA = U ook met het proces van Gramm-Schmidt bepaald kan worden. 
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Het is echter duidelijk dat zo'n proces nooit de stabiliteit van 
de Housecholder-transformaties kan halen. Bij Gramm-Schmidt trekt 
men immers van een kolom van A een lineaire combinatie van reeds 
bepaalde basisvectoren af. Neem nu eens aan dat de beschouwde ko- 
lom van A bijna afhankelijk is van de vorige. Dan zal de vector, 
die we na het aftrekken van de lineaire combinatie van reeds be- 
paalde basisvectoren overhouden, relatief sterk met ruis behept 
zijn. (ga na; zie ook (25,1)). 

Van het loodrecht staan van deze vector op de voorgaande moeten we 
ons dan geen al te grote voorstelling maken, maar dit betekent dat 
Q niet erg goed orthogonaal behoeft te zijn, in tegenstelling tot 
de situatie bij Householder transformaties. Men kan overigens het 
proces van Gramm-Schmidt modificeren, zodanig dat de orthogonali- 
teit veel meer behouden blijft. Voor détails zie Golub. 


$ 33a 
{32a.1) 


(33a.2) 


(33a.3) 


(33a.4) 


(33a.5) 


( 33a.6) 


(33a.7) 
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Stelsels niet-lineaire vergelijkingen. 


Een stelsel niet-lineaire vergclijkingen heeft de vorm 
’ fs (E1 sE) z 0 


an (Es se..sE) z 0 

waarin voor elke iz=1,..‚n Ef: R'> R. 

Schrijven we x = (E1,...,E ) en definiëren we 

fo) = CEC zemtalids sons lk p-…-»61) ), dan betekent oplossen 
van het stelsel (33a.2) het vinden van de wortel(s) van de ver- 


gelijking 


f(x) = 0 
met f : RP» RP (i.n.a. niet-lineair). 


De methoden die wortels van (33a.3) bepalen zijn globaal in 2 
groepen te verdelen: 

„- iteratieve methoden 

- minimalisatiemethoden 


Iteratief oplossen van de vergelijking f(x) = 0 betekent (vgl. 
S11e.v.) weer transformatie in (lx) = X, welke vergelijking men 
met successieve substitutie zal trachten op te lossen, 


Een afbeelding F : RP» RP neet differentiëerbaar in punt x in- 
dien een lineaire afbeelding F'(x) : RÌ > R bestaat zodat voor 
alle h: 

F(x+h) = F(x) + F'(x)(h) + elxsh) 


…_ helxsn)l 
met HIB rai 0 
oF. 
F'(x) is dus cen matrix t.wi de matrix 3E: van partiële afge= 
leiden der coördinaten van F. Men noemt dede matrix ook wel de 
Jacobiaan. 


Als a een wortel is van ò(x) = x en & is differentieerbaar in a, 
en westellenh; = x;‚-a dan geldt blijkbaar 

h 2 Xj41T0 5 bx; d-a E dlath,;)-b(a) 
d'Cadh;telh;) 


1i+1 


(33a.8) 


(33a.9) 
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met lim JeCh)t 
h0 Ih} 


zodat I$'(a)l weer zoiets als de convergentiefactor is. 


= 0 


Voor een willekeurige vergelijking f(x) = 0 zoekt mer. $ van de vorm 


bx) = x-M(x)f(x) 


waarin voor elke x M(x) een lineaire operator is zodanig dat 


lim M(x) bestaat (a de wortel van f(x) = 0). 
Xa 
Als f differentieerbaar is in a dan geldt 


Hlath)-bla) = h-Mlath)f(ath) = 
| z hM(adte, (h))(flad)+f'(a)hte,(h)) 
waarin e‚(h) een operator en e‚(h) een vector is zodat limle,(h)f = 0 
h+»0 


en bes Legt z 0. Voorts is f(a) = O0. Dus 
dlath)-dla) = [ I-M(a)f'C(a)lhtelh) 


Lel = 0. Dus $'(a) = I-M(a)f'(a). 


1 


met lim 
h»0 
Door te nemen M(x) z f'(x) “ krijgt men het meerdimensionaal' 


analogon van Newton's proces: 


(334.10) 66) =x- ft) | ex) 


(33a.11) Per iteratiestap zal men dus n° afgeleiden moeten berekenen en 


een matrix inverteren. 

Vaak is het bepalen van f'(x) al moeilijk en moet men die eerst 
benaderen alvorens te inverteren. Dit kan door in de Jacobiaan 
de differentiaalquotiënten te vervangen door differentiequotiën- 
ten en men krijgt het analogon van Koorden=-Newton (zie TA 527.u) 
Wegens de hoeveelheid werk verbonden aan het bepalen van £' 
en Era) kan het verstandig zijn dit niet bij elke iteratic- 
stap opnieuw te doen, maar slechts af en toe. Werkend met vaste 
f' Ges) krijgt men weliswaar slechts lineaire convergentie, maar 
de convergentiefactor is heel klein als X30 dicht bij de wortel 
ligt. 


(33a.12) Voor de keuze van cen startwaarde zal men zich moeten beroepen 


op zekere à priori schattingen omtrent de ligging der wortel(e), 
Indien deze ontbreken zit er niets anders op dan maar een start- 
punt te nemen en te hopen dat het iteratie proces naar een wor- 
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tel convergeert. Met een ander startpunt zal men misschien (hope- 
lijk) een andere wortel vinden etc. 

Het is duidelijk dat we op deze manier nooit weten of alle wor- 
tels bepaald zijn (indien hun aantal tenminste niet van te voren 
bekend is). Dit is ook cen heel moeilijk probleem. 


(33a.13) Bij de minimalisatiemethoden probeert men functionalen G : R- Rr” 
te vinden die kun minimum (minima) juist aannemen in de wortel (s) 
van (33a.3). 


Zulke functionalen zijn bv. 


n 
GGO = E [Eler ed! 
Ee | zr En) 
El 2 
OE 8, [£iEr oe E) 


welke nog de bijzonderheid hebben dat zij in hun (absolute) minima 
de waarde O0 hebben. 


(33a.14) Het vinden van het minimum van G geschiedt veelal met descent- 
methoden waarbij men, in de i® iterand Xi aangekomen, een rich- 
ting ri zoekt zodat G(x) < G(x;) voor vectoren Xx = X; + Ar; (met 
voldoend kleine Aà > 0). Dan zoekt men zodanig A = A, dat in de 
richting Pr; G(x) minimaal wordt, en men vindt Kia4 ° “it A;F;: 


(33a.15) Bij de methode van steepest descent neemt men 


(33a.16) r.; = = grad G(x,) 
i sG jd 

met grad G(x) zh 

E1 


SG 
lef, 
Omdat G het sterkst toeneemt in de richting van grad G zal men 
onder zekere netheidsvoorwaarden mogen stellen dat G het sterkst 


afneemt in de tegengestelde richting en (282.16) lijkt inderdaad 
heel gunstig. 


(33a.17) Toch kan het werkelijk convergentiekarakter van de methode van 
steepest descent erg tegenvallen. 
De oorzaak is doorgaans dat - grad G niet de meest gunstige rich- 
ting is waarin men het minimum moet zoeken. 
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(334.18) Bij andere methoden (bv. de geconjugeerde gradiënten methode) zal 
men vanuit Xi kleine stapjes in een aantal speciale richtingen 
ondernemen om zo direct mogelijk in het (een)minimum te geraken. 
Voor nadere -details zie bv. Ralston (ch. 8), Kowalik - Osborne. 


5 34 


(34.1) 


(34.2) 


(3u, 3) 


(34,4) 


(34,5) 


(34,6) 


(34,7) 


(34.8) 
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Eipenwaarden 


Inleiding 

Zij A cen complexe nxn = matrix. 

Een getal A < # heet een eigenwaarde van A indien een x #0 be- 
staat zodat Ax z Ax; x heet in dat geval een eigenvector bij A. 


. ie % Bd ee 
De eigenwaarden van A zijn Julst de wortels van het n graads 


polynoom Y(A;A) z= det (A = AI), het zgn. karakteristieke polynoom 


van A, 


Opgave. Toon aan cat Y(A;A) = var” AT;A) voor elke nonsinguliere 
T. 


Opgave. Zij A of B nonsingulier. 


Toon aan Y(AB;A) z YP(BA;A) . (Vla veeg  n 
von AaB sagen. En zalde gelde al, Ae Bmàr nr sa) 
VE WBAI). A en acher 


Opgave. Toon aan dat #(A;A) = W(A*;X) veracht, CWilkuason pa). 


N n 
Opgave. Toon aan dat a een wortel is van aoxì «+ ………. + An-1* * a, 


ddan als u even eigenwaarde is van 


01 
A Ss & is se " 
pe 
Ed an nd e, 0 e | 
ao ao 


Ga ook na dat (-1)7ao. VEAKA) za Bes & an- 1 + an. Afgezien 
van een multiplicatieve factor kan ieder polynoom danae optre- 
den als karakteristiek polynoom van cen matrix. 


Opmerking. Men kan Y(A;À) zien als formele veelterm en voor A ma- 
trices substitueren. 
Stelling van Cayley = Hamilton : Y(A;A) = 0, 


Indien A = A, een eigenwaarde van A is, heeft Y(A;5À) een deler 
A= A; ). Het grootste getal m zodat (A - A; jr | YCA;A) neet de 
iatneneelanes multipliciteit van A; ile. de multipliciteit 1 
is noemt men de eigenwaarde lelien tie. 


en 


(35,1) 


(35.2) 
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Als alle eigenwaarden van de matrix enkelvoudig zijn dan is er een 
basis van cigenvectoren. Dit laatste kan echter ook wel het geval 
zijn zonder het eerste (Bijv. heeft een symmetrische matrix altijd 
een basis van eigenvectoren). 


(34.9) De eigenwaarden van een matrix A zijn afhankelijk van de 


elementen van A. 

Wijzigt men de elementen van A, dan zullen de eigenwaarden 
gewoonlijk verschuiven. 

De volgende stelling leert dat de eigenwaarden 'continu' 
afhangen van de elementen van matrix A, 


elfs omalaghische) 


2 
(34.10) Stelling. Er bestaan continue’ functies Oise: Ë > Cc 


zodanig dat voor elke complexe matrix A z(a. za 


MAD = DT. Ha - 6,Can, arzo- vand) 
Bewijs.De wortels van een polynoom hangen continu van de 
coefficienten af. Meer precies: er bestaan continue 
functies &1,...; Ò: > zodanig dat voor alle 
Can s.eesa,) e &: 
n n=1 a 

Z + a1Z teeeta,? nz - Ô. elke: )), 
Ga na dat de coefficienten van W (A;A) continu van de elementen 
van A afhangen. 


Localisering van eigenwaarden. 


De volgende stelling doet cen uitspraak over de ligging der cie 
genwaarden van een matrix. 


Stelling (Gershgorin's cirkelstelling). Iedere eigenwaarde van 


A = (a. ij’ ligt in minstens cen der gesloten cirkelschijven met 


middelpunt a,. en straal £ la;;h. 
ad in j#i 
Bewijs. Zij À ew. , Xx z Ck geene dt een e‚v. bij à. 


Ui z : -; z 
it Ax z Ax volgt: (A anr) *r Ier an; xj voor alla n. 


ii 


Dan: 
|A = äepl + lx, Ee Le: ben alle n. 


Zij i zo dat |x‚| = max |x_l. Dan is x. «0. 
1 er r 1 


„ |A - ajil < e ie 


‚|. 
3 NE 


NA-112 


(35.3) Opgave. Een bewijs van (35,2) is ook als volgt te verkrijgen, 
a) kij C = les} 


ij . . Ë _ 
Toon aan dat le; ;l > LE le; | impliceert dat C nonsingulier is, 
ji i 


Vermenigvuldig hiertoe C aan de linkerzijde met zodanige diago- 
naalmatrix dat net proaukt een matrix is met louter enen op de 


hoofddiagonaal. (Kan dit?) Schrijf dit produkt als I= F 
(zekere F) en pas (19.24) toe, 


b) Pas nu a) toe op de matrix A = AI en bewijs (3542). 


(35.4) 


(35,5) 


(35.6) 


(35.7) 


(35.8) 


Opgave. Ga na dat iedere eigenwaarde van A z Ca; ;) ligt in mine 


stens een der cirkelschijven met middelpunt azj en straal 
E la; ;|- 
i#j 


De eigenwaarden van A bevinden zich dus in cirkels 


M‚(A) ={ z [iz - a;;l « AEL Ì Am Te aar atie 


Zij heten de Gershgorin _ = cirkels van A. 


ave. Elke wortel a van aak + eee + Aro4* + Ar voldoet aan 
minstens een van de volgende ongelijkheden: 


Pd 


a 
lat <1 

o a, 

lal <1 + max En 

1 bro def o 

(Aanwijzing : (34,6) en (35,b)). 


n / 
Opgave. Elke wortel a van AX tees + Arg * an voldoet aan 
minstens een van de volgende ongelijkheden: 
lal <1 


las Spr Ees! 

a+ == & . Di as | e 
ao ao! i=2 a 

(Aanwijzing : (34,6) en (35,2)) 


De resultaten uit (35.6) en (35.7) maken het mogelijk een uit- 
spraak. te geen over de positie der wortels van een polynoom 
px) = Nd + eer + ArgX ie an” 


Er zijn nog wel andere wortelschattingen bekend. We noemden reeds 








(35. 10) 
(35,11) 


(35.12) O 


(35,13) 


(35. 1u) 
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die van Newton voor polynomia met uitsluitend reële wortels. We 
noemen er nog enkele. | 

In het nu volgende zullen we steeds aannemen : ao 5 1 (geen beper= 
king). OE: 


Lemma. Laten a, t/m an positieve reêle getallen zijn met 
n 


La, < 1. Dan geldt voor elke wortel r van p(x): 


1 ia. 
el < max\ (là! 
i Ì 
Bewijs. u r 
plr) = 01 = Eek, ei k aai % _ & : | zi 
r r 1 rl 


Dan moet voor minstens één i 
la;l 
> a. 


elf 1 
d.w.z. voor minstens één i geldt 


el fi 
% 


Zo krijgt men bv. voor ai = = 2 en 


Ir| < 2 max | la, | 
i 
hetgeen vrij gemakkelijk toepasbaar is bij gebruik van een com 


puter met binair getalsysteem \Yem,2P ) met 2%” „3 ve < 29 
majorveert men dan niet al te slecht door dna +q)/1)), 


als wortelschatting 


mert Ga na wat de schatting uit. (35. 9) wordt voor 

- asl / F lajl. 5 
Geef het max. asnLietet aan. Onderscheidt daarbij tussen Ela; <1 
respe > 1, 


Opgave. Ga na dat hal g | pak tocpasbaar is met 


z /EWs, Di 


n/c. Ì 


Schrijf hiertoe e; = Vh] en beziec A Pe ° 
GC. 
J 


Wat wordt de schatting uit (35,9 ) voor deze a,. 


Het is moeilijk te zeggen welke van de gegeven schattingen de bes=- 
te is; soms is de een, soms de ander iets beter. Men kan echter 
aantonen dat (35.11) nooit meer dan cen factor 2n te groot is. 





(35,45 
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We merken op dat niet alle schattingen op dezelfde wijze reageren 
op schaalverandering. 


Bij schaalverandering gaat men naast p(x) = xr + a,xnÎ + 


bezien het polynoom met k=keer zo grote wortels. 
Ga na dat dit het polynoom 
n n=1 Kri 


q) = Xx + k ax + ves + 


n 
1 xXx + k a. 


An-1 n 


is. 

Als de schatting voor de wortels van q(x) kx zo groot is als de 
overeenkomstige schatting voor de wortels van p(x) noemt men de 
wortelschatting homogeen, anders inhomogeen. 


(35. 16) Opgave. 


(35.17) 


(35.18) 


a) (35.44) en (3543 ) zijn homogene schattingen. 

b) (35.42), (35.6) en (35.7) zijn inhomogene schattingen. 

Ga voor verschillende polynomia met bekende wortels (bijv. 1,1,1 
en 106,10%,10°) eens na wat er uitkomt. 


We formuleren nog een tweede stelling, die een 
interessante eigenschap van de Gershgorin=cirkels 
geeft. 


Stelling (Gershgorin's ze cirkelstelling ). Indien k 
Gershgorin-cirkels van A geïsoleerd liggen van de 


overige, dan bevatten deze k cirkels tezamen precies 

k eigenwaarden van A, elke eigenwaarde geteld naar 

zijn multipliciteit. 

Bewijs: 

Schrijf A = D + B, D een diagonaalmatrix met djj=ajj dzlseeeon). 
Zij At) = D + tB voor t € [0,1 1, dus Alo)=D en 


AC1)z A. Zij C‚(H) = {z | Izea;; | S t Za! a;s! }li=lsees on) 


Dan zijn de c‚ Ct) de Gershgorin=cirkels van Alt). 
Wegens (34.10) bestaan continue functies Xjseee kn: R > € 
zodat XC), seo Xr (HD de eigenwaarden van Alt) zijn. 


Na eventueel vernummeren der functies X; geldt voor 
alle i: X;(o) = Ajj* 
Stel nu dat k Gersh gorin=cirkels van A geïsoleerd liggen 
van de Overige, Het is geen restrictie te onderstellen 
dat dit C,CID jee CCI) zijn. Evident liggen dan voor 
alle t& [sil CH), Gt) geïsoleerd liggen van 

de overige c‚ (t). 


(35,19) 


(35.20) 


(35.21) 
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Als nu b.v. XC) 5 C41(1)s dan moet wegens X, O)za,, 


een t bestaan zodat Xx, Ct) in geen 
enkele C‚ (1) ligt, dus zeker in geen 
enkele C‚ (t). (zie b.v. nevenstaande 
figuur en (35,19)). Dit is een tegen- 
spraak! Etc. 

Men ziet in dat XD ses oC) 
noodzakelijk in de vereniging der 
C4(1),...; G, (1) moeten liggen, 





In de bovenstaande figuur is het geval geschetst dat van 
reele 3 x 3 A twee Gershgorin=cirkels geïsoleerd liggen 
van de overige. 

De aangegeven cirkels moet men loodrecht op het vlak van 


tekening denken. 


Beziet men de Gershgorin=cirkels met middelpunt ai 


t van 0 naar 1 loopt, dan vormen zich kegels, waarvan de 


als 


projecties zijn aangegeven. 
De functie x, zoals getekend kan nooit het verloop van 
een eigenwaarde zijn. Die moet te allen tijde binnen een 


Gershgorin=cirkel Liggen, hetgeen betekent dat xt) 
op elk moment binnen een kegel moet liggen. 


Opgave. Ga na hoe het argument uit (35,19) in een 
willekeurige situatie te voeren is. 


Opgave. Zij A= D + €B, D een diagonaalmatrix met 


des 8 Ass 
11 11 


Stel voor alle iz 2,...‚n : ai Ea az zij é de: 
Toon aan dat voor alle e met lel < Tier binnen de 
Gershgorin=cirkel rond as1 precies één enkelvoudige eigen=- 
waarde van A ligt. 


5 36 


(36,1) 


(36.2) 


(36,3) 





Perturbatie van eigenwaarden. 


Verstoringen in de matrix A beïnvloeden de ligging der eigenwaar- 
den en eigenvectoren. 

we zullen dit effect onderzoeken voor het geval dat A een diago- 
naliseerbare matrix is. 


„A e 
Laat ri AT = D, met T nonsingulier, D | , ) een 
“An 
diagonaalmatrix waarvan de diagonaalelementen dus juist de eigen- 
waarden van A zijn. 
we bezien eerst het effect van perturbaties € B, 


Ee 
Zij T-” BT = (8,5). 


Stelling, Zij A; enkelvoudige eigenwaarde van een diagonaliseer- 
bare matrix A, 
Dan heeft A + eB voor voldoend kleine e precies éân, eveneens enkelvou- 
dige eigenwaarde H‚ waarvoor geldt: 
Ius = A, -e Bijl Sk. lel? Rn 18,5! 
met zekere constante «,die onafhankelyk van £ is, 


Bewijs. 
De En van A + eB zijn gelijk aan die van 
T KS + EB)T = D + ET 5 BT, en dus ook gelijk aan die van 

1D + eT”Î BT)P = D + e P“Âr"1 g T P voor willekeurige nonsin= 
Aer diagonaalmatrix P. 
Kies nu in P het i° diagonaalelement gelijk aan (ker Î, en de overige 
diagonaalelementen 1. 


pipi Br p ontstaat uit T 


il 


“1 pr = (B, 3’ door de i® kolom te ver- 


menigvuldigen met (ke) en de i® ‘erij cn vermenigvuldigen met ke. 
Het element op de plaats (i,i) verandert dus niet. 

Bezie ‘le Gershrorin = cirkels van D + eP“ÎT”Î BT P. 

We weten uit de enkelvoudigheid dat A, verschillend is van alle o= 
verige elementen op de heetaätanshaal van D. Ofwel: de Gershgorin- 


cirkel M; (D) (straal = 0) ligt geïsoleerd van de overige Gershgorin- 


cirkels. ekere 

Maar dan zal voor e met |el rg Gershgorincirkel van 

D + eP ir ÎBrTP met middelpunt A, + eB‚; en straal 

Kle}? 5 |B, 3! ook nog geïsoleerd ann van de overige, mits « 


voldoedäd an en T voldoende klein. 


(36.4) 


(36,5) 


(36.6) 


36,7) 


(36.8) 
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Wegens (35.6) zal deze cirkel juist één, enkelvoudige eigenwaarde 
van D + ep”dr”ÎB T P, dus van A + €B bezitten. 

Uit (36.3) ziet men dat een enkelvoudige eigenwaarde Aj van een 
diagonaliseerbare matrix A in eerste benadering verstoord wordt 
met een bedrag €B;;° 


Opgave. Toon aan dat men in het bovenstaande bewijs kan nemen: 
| + is 
Elenjl + 18,5 D 
1: È 
eg Ary Acngorolis eordoara, walnicas 
Gaat men het effect van perturbaties €eB op een niet=enkelvoudige 
zeg m-voudige eigenwaarde A na“;dan vindt men analoog aan het 
bewijs van (36.3) m Gershgorincirkels die voor voldoend kleine 
e geïsoleerd liggen van de overige. 


Dit geeft niet veel meer inform:tiec dan dat er m eigenwaardenu van 
A + €B bestaan zodat 

Iu - Al Se. IT BTÍ (Ga na). 
De volgende stelling leert dat algemener geldt: 


d 


Stelling. (Bauer = Fíike), Zij A even diagonaliseerbare matrix, 
T° AT =D. 
Zij u een cigenwaarde van A + B, 


Dan geldt: min [u - A; | s Iri BTH, voor alle pelt‚ael- 
i 
Bewijs. We mogen aannemen dat min| u - A; | * 0, dew.z. dat u geen 


eigenwaarde van A is (anders wa er niets te bewijzen), 
TTÎ(A - LIJT = D = uI is dus nonsingulier, 
TA + B = UIJT = D= uI+ Tí B T is singulier zodat ook 
D- ui eDp-ureTÄ BT) =I+De ui iT ÄpT singulier is. 
Wegens (19.24) moet dan H(D - ui) ir” BTI_ > 1 voor 1 < p S ee, 
Ga na dat I(D - Dhr erh, < nT Tr” Apr D* 

i 


Voor symmetrische A volgt met p = 2 dat binnen een afstand MBI, 

van een eigenwaarde u van de gepertvrbeerde matrix A + B minstens 
één cigenwaarde van A te vinden is. Door deze stelling wordt even- 
wel niet uitgesloten dat alle eigenwaarden van de geperturbeerde ma- 
trix in de buurt van Één bepaalde e.w.van A liggen. 


pnt ne 


(36.9) 


(36,10) 


(3511) 


(36.12) 


(36.13) 


(36.14) 


(36.15) 
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Als A en B beide symmetrisch zijn kan dit laatste niet optreden 


en is een veel preciesere relatering van de eigenwaarden mogelijk 
ie echter 36 13) 2 


Stelling. Zij A symmetrisch met eigenwaarden Ai < Aa &... < An: 
Zij A + B symmetrisch met eigenwaarden ui S uas... < Ho 

Dan: Lu; hd A; | % 1BĲ, s Ì = 1, eee Ne 

(Voor een bewijs van (36.9) wordt de geïnteresseerde lezer verwe- 


zen naar Wilkinson - The Algebrafc Eigen value Problem, pg 99-102.) 


Opgave. Zij A = (ass) een symmetrische matrix waarvan de coëffi- 
ciënten niet alle computergetallen zijn in de zin van $ 4 , De k° 
eigenwaarde van A en de k° « eigenwaarde van de computer=represen- 
tatie van A verschillen hoogstens n. VAL . E‚, waarin 


bal = max |a, 


kad ij 


Opgave. Zij A = (a;;) een symmetrische matrix met eigenwaarden 
A1 Sha So. SA 
Zonder beperking mogen we aannemen: art S a32 & .. < am" 


Dan geldt: NA; - a;;l < max Er la;:l. 
De 1 JA 77 
(vgl. (35.2)). 


Zij A symmetrisch, B een symmetrische perturbatie van A. 

Zo'n perturbatie kan bv. het gevolg zijn van afrondfouten in de 
coëfficiënten van A of optreden als de coëfficiënten empirisch 
bv. door meting bepaald worden. 


Voor het verschil ÀÀ, = - Ar (u en A als in (36.9) ) geldt: 


k 5 Ek 


ze zi Ed 
|A) & HBI, = oC) &n. KBh, 


zodat wegens p(A) = lAÁ, > max |a. 


ijl = MAI, geldt: 


ed 


es 


| 1BI, 
Sn. TAN * 
) mm 


ä 


De verstoring in de k° eigenwaarde relatief t.o.v. de absoluut 
grootste is dus hoogstens n=-keer zo groot als de maximale ver- 
storing in de matrixelementen relatief t.o.v. het absoluut groot= 
ste matrixelement.Dit sluit natuurlijk niet uit dat 

Lan, mn 
Dl best groot kan zijn. Men lette bv, op de eigenwaarden die 


in abs. waarde klein zijn t.o.v. p(A). 
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(36.15) Voor de verstoring in de k° cigenwaarde tengevolge van afrondfou- 


ten in de coëfficiënten van A geldt dus (zie ook (36.10)): 

hea Oe 

PA) * Ne E 
zodat men mag stellen dat ven symmetrische matrix zijn eigenwaar- 
den op een uitstekende wijze representeert. 


(36.17) Bij niet-symmetrische matrices kan de situatie vecl ongunstiger 


zijn. 


(36.18) Voorbeeld: 


20 2Q 
De À 
A= | Bes 
\__@ °° 20 
“4 4 
0 
Bef £ (20 x 20) 
onm e ed X 


De eigenwaarden van A zijn 1 < 2 < ... < 20. 
Y(A + EB3A) = “Ml (ied) - 20ÎPe = W(AIN) - 201% 
iz1 
Omdat op [1,20] : FYA; DD | < 20: (ruime bovengrens) zal voor 
€ > 20e, “ U,‚5 10”? het polynoom Y(A + eB;A) nog slechts 2 reële 


wortels hebben, de overige zijn complex! Voor es}, geor ondeha. vore naded, 
(36. 19) Stalles Zj A Aagomalisanbaan . Zi B wllelaur Dan bestat er tw 
Namrmararnag Ve da. Lgemumsndon A va A en po ve. AFB zodat [As -uol 
(em-4) IUT-!BT Ien , T de deagemaliaurwnde Dn rmdhe vaa. A. 
Bows Tola de A và Aat Crvplear vlak … Vd elke Frera 
den Ac olà mm abarlade woeards Ameen dae, 2 IT BTI va, learn A 
baage. Benchruas en, Hue van de Zo onbdauu Eje, 
grnap De bbrshpusoirdala Kvvean Doggen, geiseteund vaa, De 
te hendel ae ada A's, Nu 
Inde, p h % ; 
henk Mmmm Wer geochilde Mumruig dan 


"Ee 


CA, 4) 


637, 2 


(473) 


(37.4) 


(37 A) 


(37.6) 


(37,7) 
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Perturbatis van eigenvectoren. 
ne neden nt Et ed enten 
Zij A weer een diagonaliscerbarc matrix en stel 


{u Ì 
-1 j 


ae on 
 - Sa 


Zij A, een enkelvoudige eigenwaarde van A. 


T 


9 
a 


n / 
Het is geen beperking i = 1 te neinen. 


Uit (36.3) volgt dat onder alle eigenwaarden van A + €B met 

lel < tT, tT > 0 voldoende klein, er precies één, zeg ule), is die 
het dichtst bij A: ligt. 

Deze u(e) is bovendien enkelvoudig zodat de bijbehorende eigen- 
ruimte evenals die van A, 1=dimensionaal is. 


Zij vle) een cigenvector van A + eB bij eigenwaarde u(e), v een 
eigenvector van A bij eigenwaarde Ai (vle) en v zijn op scalaire 
factoren na uniek bepaald). 


We willen nagaan hoe goed de richting van v(e) de richting van 
v benadert. 


Gebruikmakend van de SADE ASC BP RESEREAN gaan we natuurlijk kijken 
naar z = T dy en z(e) = ri vle) die eigenvectoren van T° * AT = D 


resp. TÎa + €B)T = D + e(B,;5) AAN. 


Lemma. Zij C een singuliere n xn matrix die gepartitioneerd is 
als (P fe) ‚ P een k x k matrix, S een (n=k) Xx (n-k) matrix. 

zi Ee Ogseees Xn )f een vector waarvoor Cx = 0. 

Partitioneer x als , waarin u staat voor de eerste k coördinaten, 
v voor de overige (n=k) coördinaten van x. 

Zij S nonsingulier. 

Dan geldt v = =S°ÌÎ Ru. 


Bewijs, Ru + Sv = 0. 


We passen dit lemma toe voor k = 1 op de matrix 
Ca B 4 e(B,5) e HE) La 


(37.8) 


(37.9) 


(37.10) 


(37.11) 


(37,12) 


(37,13) 


(37,14) 


(37,15) 
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Lemma. Voor voldoend kleine e is de (n=1) Xx (n-1) rechtsonder- 


matrix S van C = D + e(B;;) — u(e) I nonsingulier. 


Bewijs. Definiëer Kk = lr Ss As! 


Zij T1 > 0 zodanig dat voor alle e met Jel < |ul: 


min fule) — Al > jk. 
it 
2 max E |B; +! 2, AT “sThe 

im j=t *J 

Ga na dat voor del < min(t,T: ,T2 ) de matrix S voldoet aan de con 


dities van (35,3a) 


Derhalve kunnen we voor de eigenvector z(e) van D + e(B,5) stellen 
dat 


u\ 
zie) = (8) 
-1 

met v z =S R u. 


Hierin is u in feite de eerste coördinaat van z(e) die we gevoeg- 
lijk gelijk aan 1 mogen nemen. 
Dan is 


ze) ( 


met 


(eBa: vts” 


dat voor voldoend kleine € de vector z(e) continu afhangt 
dat in het bijzonder z(e) steeds meer de richting van z 


Men ziet 
van € en 
krijgt. 

D' + eB' 


Schrijf S = met 


- 


‚Aa-ufe) 
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B "Tr" Kgaî 
Ei | 
! , 
| 
Lj 
n2” Ne MS Brin/ 


Dan is wedens (20.42) ev w efens pede A4 HOE) 


Ü 
i ene 
nn eN 
Nn Aa ple vote) Tha, + Ole) 
S= ge An A, 
An mle) ki 
zodat fn 
ve & [Ze \ zoen 
Éna 
As =An 


(37.21) Voor voldoend kleine e mogen we derhalve stellen dat in eerste 
orde benadering: 


€g €8 
(37.22) Ee anr ree he 
1 A3 adds \ 
Terugtransformeren met T geeft dan de 


ve). 


2E SCH, 


eerste orde benadering voor 


r 


(37.23) 


(37.24) 


EMED 


(37.26) 


(37.27) 


(37,28) 
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Laten we nu aannemen dat T unitair is. 
Ga na dat dit juist optreedt in het geval A een normale matrix is. 
(vgl. (19.10) en (19.12) ). 


Omdat T en dus Tr! inprodukten behoudt, is de hoek $(e) tussen 
ve) en v gelijk aan die tussen «íe) en z zodat: 


n [8 A Es È n 3 
tg dle) > |lel.{ £ he: EG <4 [8,,1° . fel 
j=2 Dag-A,l 2} 


We hebben nu: 


Stelling. Zij A; een enkelvoudige wvigenwaarde van een normale 
matrix A, Kk = minlA,=A. | 

El 1 J 
Zij € voldoende klein. 
Onder alle eigenwaarden van A + eB is er precies één (u, ‚(e) ) die 
het dichtst bij A; Eek 
Deze u; Ce) is enkelvoudig en voor de hoek $(e) tussen een eigen- 
vector bij u‚(€) en een cigenvector bij A, geldt: 


|ote| & Eik le | 


De perturbatie van cigenvectoren bij een meervoudige eigenwaarde 
is een aanmerkelijk moeilijker aangelegenheid. 

Een eigenvector is dan nict meer op een scalaire factor na uniek 
bepaald. 


‘Men kan hoogstens zeggen dat de eigenruimte bij de gepertubeerde 


eigenwaarde steeds dichter "nadert" tot de oorspronkelijke eigen- 
ruinte indien e *0, 


Voor symmetrische perturbaties A + eB van symmetrische matrices A 
kunnen we dit nog nader analyseren, 

De volgende stelling is ook toepasbaar op het geval dat een aantai 
eigenwaarden van A dicht bijeen liggen of samenvallen. 


Stelling. Laten A en B symmetrische matrices zijn. Stel A heeft 
eigenwaarden Age ve » Àn (in een of andere volgorde), 

Zij V de eigenruimte opgespannen door de eigenvectoren bij 

Agere A EN zij V het orthoplement van V. 

Zij A de geperturbeerde waarde van een der Ajs Bee 3 Ar zodat dus 
|a-Ajl S lel WBEa (zekere iz1,...k) on zij x een bijbehorende 
eigenvector van A + eB, 


Schrijf Xx z Xx + a ú € Vv, ne € vt 


(37.29) 
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Bij es min [Az =2elBl, > 0. 
Jk, 1% 
han is let, < Lel IB, belt. 


Bewijs. De stelling is basis=onafhankelijk en we mogen dus ver- 
onderstellen dat A nevenstaande gedaante /À4 | \ 
heeft. Partitioneer A + eB -Alevenzo: Á gn 

EE Q) . Dan zijn alle e.w. van S in abs. Al 

waarde > ps dusts"il, < À. voorts geldt | — — je == 


o k+1 
x x Ee 
k ka), -s"tr[:° | | hk 
5, ) XL 

waaruit wegens ÍRI, S lef IBE, het beweerde volgt. 


Bijgevolg, als de eigenwaarden Agsveed van A weinig van elkaar 
verschillen of eventueel samenvallen, maar behoorlijk gescheiden 
liggen van de overige eigenwaarden, dan heeft de bijbehorende 
eigenruimte van A + e€B een basis waarvan alle vectoren hoogstens 
een hoek » 5 IB, met de eigenruimte bij Agseees A van A maken. 
Ben mog vond brharfur Maulbaat ralf mtan vi het Kerr foaaie (moor moedpfa) 
arbchil van C. Darts en WM Kohân 1 Tt zetehdn U eeewectrro ba 
perZunbaien (MT , LAAN. Num Brad. FUR) 1-teb: 


(3 30) teluig. Leds PoB Epnmebrischt mulier em. GW Cepeda 
behorend 


dy le Oegenrvnan dan Agens Az Ute, (oa ta of ader 
Tiemarkrag ) Katt. As, …Âk hi U. ordt (2, ) Lg n le Aloys en Äm 
er buil. riana Pei ueklef ikan ded 
bize rn uwer behoud de] ole eigsmras wl dl, HreB Od we U Zen 
Ven ( BG) Es A, Ig Krt. Z Lb mitemal heek Zouden A beun 
Fn W on haer fragiel op V Dan wo Sl) < FIB 
Bumgs (eeu voudleger da. by Daver- Ka han) Asa Thet Loor 
dlbktunige % Art u la vga afs AAA 5 dare ple 
Bfurhrg aa B Want dak Akelig AoA Le Sepmank [-ÀF] 
Ofafrarr nen, heli Aers s … An, Ak Hema het deralk aa Ht anders 
Amadoeg Mk Kagan dass Le slade Via da belang, UManù Pap 
n= <1IBĲ, _ 0 


e- | ad) 


he Me add 


Alen st. Am 
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A beet VE ps zeekrelf nf, e. AAA pelt Wia (+ pen) 
(zä (20.6), zeelar hx afpaldeng Pp ed. | 

AseBhertolt W ú zoekzalf RP, On Clsrap geldt MBE B) IS (a) 
(ze (2e.6) e(36.5)). 

Äj Mau u, w eem Veche Jaar, uevd, web, J|ull= lwil=1, 
Wai (kw) Crijnoduct) mattaal bs. Bt ten VE VL, Ivil=1, 
zorlat Av=cu, CPO, ondud eZtp+n Ma geldt Ertunijds 

| (Av‚w) — (@+£8)v, w)| = €I(Bv‚w)l sE NBlj= 7 

Uv) -(A+EB)v,w) = (Av, w) — (Vv, Are B)w)z 

(ftp an) luw) (Yar) (uw) = P (WoW) 
zedot (Uw) < B , 


Zj ma xeW x wilbikuurig  bhaf renem, me Vane VL 


Dan cs PA(L(M Ae en(t(rma) E(wwds Yp. 
Be, SB) < GIG. 
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3 38 Methoden voor de bepaling van eigenwaarden en eigenvectoren; 
inleiding. 


Er zijn 3 typen van methoden te onderscheiden: 


=ilethoden waarbij het karakteristieke polynoom v(A;A) expliciet 
bepaald wordt. (De vigenwaarden vindt men dan als wortels van 
dit polynoom). 

“Methoden dic berusten op rechtstreekse bepaling van de Bulpun- 
ten van de functie y(A;a) = det(A-AI), maar daartoe A eerst 
transformeren tot een gedaante waarvan gemakkelijker de determi- 
nant te berekenen is. (Methoden van Householder, Givens) 

“Methoden die berusten op een rechtstreekse iteratie met matrices 
(machtsmethode, Jacobi's methode, LU en QR = algorithmen) 


5 39 Methoden die het karakteristieke polynoom expliciet bepalen, 


(39,1) Een uitgebreide behandeling van deze methoden vindt men in 
Faddeev & Faddeeva (Ch IV). 
Zie ook bv. P.A. White Jsiam 6 (1958) 393-437, 


(39.2) De methode van Krylov maakt gebruik van het feit dat Y(A;A) = 0 
(zie 64.7) ). 


n Nn 
Stellen we Y(AsÀ) z E ar A met a, = 1 dan geldt voor wille- 
keurige x: 0 p 
V(AZA)X = (2 aA) x ss Û. 
8) 
Genereer nu de rij {x,} iz0,...‚n Met 
Xo = X 
ek 
dan geldt: Dn 
La s 4 
0 k “k 


en dit geeft in coördinaten een stelsel van n lineaire verge- 


lijkingen voor de coëfficienten las bsao,... ‚ne1. 


Ga na dat het bepalen der coëfficiënten ons op deze wijze 


ed zn’ AV kost. 


(39,3) Opgave. Ga na dat het stelsel singulier is als A twee Jordan= 
kastjes heeft met gelijke eigenwaarden. 


(39,4) 


(39.5) 


(39.6) 


(39,7) 


(39,8) 
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Andere rekensehemna's zijn afkomstig van Leverrier, Danilevski 


a Cen hieruit 


kan men de coëfficiënten van het eigenwaardenpolynoom zonder 


v.d, en berusten er bv. op dat spoor (AX) z= EA 


veel moeite bepalen). 


Vanzelfsprekend maakt men bij het bepalen van de coëfficiënten 

van het eigenwaardenpolynoom afrondfouten, soms zelfs heel ernstige 
(vgl. (39.3) ). Maar zelfs als men de coëfficiënten heel nauw- 
keurig te pakken krijgt, kunnen de: wortels nog ernstig verstoord 
zijn. 


Als afschrikwekkend voorbeeld van de invloed van kleine pertur= 
baties op de ligging der wortels moge dienen het karakteristiek 
polynoom van de mapix A uit (36,18): 


YOD == MA (Asi) z ARO 4 arol? +, + Be « 
iz1 
Aangezien de wortels van een polynoom analytisch van de coëffi- 


ciënten afhangen geldt voor een (enkelvoudige) wortel p van Y: 


k k 
d = (-1)P*! 
dps Ar Gas CDP end Aa, 
Voor p = 16 en k= 15 : Ap * 4.10° Aa: 


Stel dat a, s (*1010) in de rekenautonaat waarmee men de oplossing 
wenst te bepalen niet exact voorstelbaar is, 
De afrondfout a, 8 (vgl. (4.9) ) veroorzaakt een verschuiving 
van de wortel p = 16 met 

Ap =u,10', a,s 8 u.101* E, 
In het geval &£= Ez zg = 310712 (vgl. (4,10) ) is dit onge- 
veer 200, een bedrug dat gezien de separatie van de wortels in 
de oorspronkelijke veelterm veel te groot is om de gepleegde cer- 
ste orde benadering in (39.6) tc rechtvaardigen. 


14 


Het blijkt echter wel duidelijk dat bij de polynomen met een der= 
gelijke verdeling van nulpunten alleen al de afronding die nodig 
is om het polynoom in een rekenautomaat met een rekenprecisie van 
12 decimalen te representeren de wortels onherkenbaar verminkt. 
De bepaling der eigenwaarden van een matrix met eigenwaarden in 
de buurt van 1, 2, 3, .….…. , 20 d.m.v. de oplossing van het expli- 
ciete karakteristieke polynoom is dus op een dergelijke machine 
niet mogelijk. 

We wijzen er nog eens met nadruk op dat dit in laatste instantie 
niet veroorzaakt wordt door de afrondfouten tijdens de bepaling 
der coëfficiënten van. KA;A) (die maken de zaak alleen nog maar 
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erger) maar door de onmogelijkheid de coëfficiënten adequaat in 


de machine te representeren. 


(39.9) Gelukkig zijn er heel wat methoden die volgens andere principes 
te werk gaan. 


(UO, 1) 


(UO. 2) 


(LO,U) 


(HO, 5) 


(40,6) 


(LO, 7) 


(H0,8) 


(10,9) 


NA-128 
De machtsmethode, 


A1 hd 
Zij A diagonaliseerbaar, T”ÎAT = Dl 
Neem aan dat [Al > |A, > …. > Danl. 
Zij Xo € X een (in principe) wil tekeurijs vector, 


X heeft een basie Vjseses Vn van eigenvectoren van A (vj eigen- 
vector bij A;), dus we kunnen schrijven: 
Xg 5 azv4 + avs P eee + Ann = a1v4 + tg: 


Wij bezien de rij geïtereerden Cy} "met 
in 


Xe S Axyng (k > 1), 


Dan: 


k, 
nn 


Xe = ak X9 = arAivs + eee + andnv 
AK (aqva + Tie) 
1 11 k 


met 


LL 


Az.k Ank 
Tk 2275) Vg + cecee + an Vn 


À . 
Omdat || < 1 voor i #1 geldt lim 1 = 0 en zal Tt, > 0 en dus zal 
kj k k 


Xj Voor k > ee steeds meer de richting van Va de "eerste" eigen 
vector aannemen (mits a, * 0), 
Er geldt: 


…_ Mja 4l 
lim rr = Dal 

…, (Xke1, Xk) _ 
lim Kes A4 


In de praktijk zal men Xk door een geschikte factor delen om 
te voorkomen dat de coördinaten van X onbepaald groot of 
klein worden. 

Het is gebruikelijk de geïtereerden te normeren op 1 

(bv. in de 2- of s=norm), 


Xx 
Ook y‚ = TT neemt steeds meer de richting van vj. aan 


… lAyyl _ 
en lim Tyr * ha, | 


… (ÀYk‚yk) _ 
TP Ge HM 


(HO,11) 


(40,12) 


(40.13) 


(40,14) 
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Uit (40.4) en (40.5) ziet men dat de convergentiesnelheid 
van het proces bepaald wordt door de verhouding Ee; indien 
|A,l = |Ayl kan de methode onder omstandigheden erg langzaam 
convergeren. 
Voor manieren om de machtsmethode te versnellen zie bv, 
Wilkinson (Ch 9), 


Heeft men A, en een bijbehorende eigenvector vj gevonden 
dan zou men met een startvector xy Ll vy in principe As 
(mits [Al > |Ajl voor i > 3) kunnen vinden, 

Iteratie met zo een x, zal echter t.g.v. afrondfouten spoe= 
dig een vector opleveren die toch een component in de 


richting van v4 heeft. 


Om toch de overige eigenwaarden (en eigenvectoren) van A te 
vinden zijn diverse deflatie-technieken ontwikkeld, 

Daarbij wordt na bepaling van Ay en vy op speciale wijze 

uit A een matrix B geconstrueerd die eigenwaarden 0, Asse» 
An en dezelfde eigenvectoren vjses«»Vn heeft. Etc. 

Men moet wel bedenken dat steeds gedeflateerd wordt met be= 
naderde eigenwaarden en -vectoren, hetgeen met name de lig- 
ging der nog te bepalen eigenvectoren aanmerkelijk kan be- 
invloeden.(Het effect op de nog te bepalen eigenwaarden 
valt meestal wel mee). 


De grote hoeveelheid werk en de niet erg grote nauwkeurig- 
heid zijn er de oorzaak van dat de machtsmethode thans niet 
veel meer gebruikt wordt, vooral omdat men nu over betere 
algorithmen beschikt. 


Een aan de machtsmethode verwant proces is dat der inverse= 
of Wielandt-iteratie. 

Stel men heeft op een of andere wijze een benadering \s van 
A; gevonden. … 

Als men aanneemt dat A; dichter bij A; dan bij enige andere 
eigenwaarde ligt dan is Aj. A; de absoluut kleinste eigen- 
waarde van A = AGT zodat (A; = At de absoluut grootste 
eigenwaarde van (A = A1 is. 


Men kan nu in principe de machtsmethode gaan toepassen op 
(A = Rele, maar om het onvoordelig inverteren te omzeilen 
(vgl(22.21)) gaat men als volgt te werk: 


5 ui 


(41,1) 


(41.2) 


(u1.b) 
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Kies een ao _ 
bepaal xj als oplossing van (A - AzI)xy = Xg 
bepaal xj als oplossing van (A = A12 = X4 
ate. 
Act proces blijkt vaak bijzonder snel te convergeren en men 
kan aldus nvy cen verbetering van de oorspronkelijke Ks verkrijgen. 


Het gedrag van matrix » vector =- iteratie. 


Informatie omtrent het gedrag van de rij vectoren Akx(k z 0,1,2,..) 
is van belang voor tal van (iteratieve) processen (zie bv. 528 
539 en 840). 


Het is om deze reden dat wij er nader aandacht aan schenken, 


Zij Xp * X 
XT Ak 2 1) 
Laat B z ri AT de Jordannormaalvorm van A zijn (vgl. (19,16) ), 
en definieer: 
Yn = Tis 
Y = B yi * T Xy k> 1 
We mogen aannemen dat de diagonaalelementen Ag EEN An van B 
voldoen aan LA, | > hAl B avan oe 


Beschouw eerst de situatie dat B uit slechts één Jordan-kastje 


bestaat: 
À 1 
/ N " N 6- 
N. 4 \ 
B el 8 (n xn) 

\ Ne 
\ Es ° 
A \ 1 
6 | ki 

z À I+N 


Toon aan dat NP = 0 en dat 


ak A kel k \ak=n+1 
‚An ek ee Hei) \ 


1 Hi NN 
hor 


NA bn Ì 


(11.6) 


(41.8) 
(41.9) 
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Zi nu Yp * (Es de ET en zij 1 de hoogste index waarvoor 


Aj #0 


Dan heeft Yx 5 gk Yg Als verste coördinaat: 
ZN \ 
k {k\ ‚kel - k\ ‚kel+1 
ha ad Ls SN el 
Voor k * ee is alleen de laatste term hiervan belangrijk en men 
kan nog zeggen dat deze gelijk is aan 


1-1 
kriel E, cavol£)o. 


=1): 
Evenzo gedraagt zich voor k * ee de tweede coördinaat van Yi 


k=-1+2 1 
À ejcneo(f)p 
en is dus klcin t.o.v. de eerste coördinaat. 


Idem voor de overige coördinaten. 
Ga na dat men nu kan schrijven : 


k 


als 


fe} 
_ vl=1k E) 
Ye = k ir « off 
\0 
waarin I de eenheidsmatrix is en Te) een matrix voorstelt. Voorts 
Ë 5 
es Zyaitt, (1-1) 


Laat nu B uit meerdere Jordankastjes bestaan. Stel dat het is 
Jordankastje de afmeting n, xn; en als diagonaal element A; 
heeft. 

Partitioneer Yo OP overeenkomstige wijze als B: 


ie T 
Yo had (E44 ee San,” Eier San,’ sne) . 


Zij 1, de hoogste index waarvoor Ein, * 0 is, en zij 
51, f 


_ | 
cj zn C1j=1)t (vgl. (41,8)) 
1 


Dan geldt: 


(461,10) Vk =k 


(61.11) 


(L1.12) 


(L1.,13) 


(u1.1u) 


(H1.15) 


(41.16) 


(u1.17) 


(H1.18) 
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| 


A 
mad 


_ a 


ks 1 


ALC + of)» 


5 
amd 
men, 
a 
re 
eon 


I 

! 

L 

| 

| 

| 

|. 
en 
> 

a | 
| 
' 


. eee 





/ 


otret me dat de „uPste ) Jordankastjus allemaal dezelfde eigen- 
waarden hebben, dus Nr = Az Z eee 5 Aj en A | > Ds | 
en laat ec, t/m e: niet allemaal O zijn. Dan spelen de overige 
vakjes praktisch niet meer mee. Er geldt: 

l,=1 l,=1l 
: A5 (II+ oi ces; Bins 2 


Buran GE 
Ag EE & of Yke1 


Et 
1 
Y 5 k Cs EE les: 


We bezien wit dit betekent voor de oorspronkelijke matrix A. 
Via (41.3) ziet men: 
1 1 


á n ae - 
Ke Ty tT e ol) TÌ ns, = AIT + (EN ens: 
en derhalve: 


bx, „5 [| EE 


Voor k * 99 verandert de "lengte" van ven iterand bij elke stap 
met een factor die ongeveer gelijk is aan de spectraalstraal, 


Opgave. Ga na dat onder de gemaakte aannamen Y en dus xj, Voor 
k eo tot een zekere limietstand nadert, 


Opgave. Ga na dat (UL) t/m (41.17) ook nog gelden als de eigenwaarden in 
eerste j Jordankastjes verschillend zijn maar wel dezelfde modulus 
hebben dew. 2. DD, |=... = |A: | en |JAsl > Dj |, mits onder 

de ys ves 1; precies één erootste is. 
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(h1.19) Opnerking. Tengevolge van afrondfouten in het iteratieve proces 


zullen de 1, vrijwel steeds hun maximale waarde hebben d.w.z. 


ls S ties 
1 Ì 


(t1.20) Opgave. Ga na dat voor elke matrix en elke beginvector Xp geldt 


Meh = Dsl tx, voor zekere eigenwaarde À (afh. van Xg)- 


NA-134 


Lass 


Methode van Jacobi 


(42,1) Evenals de machtsmethode behoort de methode van Jacobi tot de 
derde categorie van $38. Á 


(H2.2) Zij A = (a; 5) een reële, symmetrische matrix. 
A kan door een orthogonale transformatie op diagonaal vorm ge=- 
bracht worden d.w.z. er is een orthogonale O zodat 
otaoO=De 0 © ‚ met D een diagonaal matrix waarvan de 
6 ‘vAn 
diagonaalelementen juist de (reële) eigenwaarden van A zijn. 


(412,3) De door Jacobi in 1846 aangegeven methode beoogt 0 successief op 
te bouwen als produkt van rotaties in geschikte vlakken. 


(42,4) Er wordt gebruik gemaakt van de volgende eigenschap die we als 
opgave formuleren. 


(42.5) Opgave. Voor een willekeurige matrix C en unitaire matrices U 
en V geldt tucvl, z ich. 
Bemerk hiertoe dat een kolom van UC dezelfde 2=-norm heeft als de 
overeenkomstige kolom van C en dat een rij van CV dezelfde 2-norm 
heeft als de overeenkomstige rij van C. 


(42.6) Zij nu voor p <q: 


he eG n 
De “eur--s DP 
Pa Hd “1 1 e 
GQ … ee Ce K mmm | 
8 1 î 1 
& e 
p= q= 
met c? + s% = 1 


(42,7) Dan is Ora orthogonaal en in feite een draaiing in het door de 


basisvectoren ep en eq opgespannen vlak. 
_1 T 
: O0 z 0. 
DU Spa Ara 
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(42.8) Zij a een absoluut grootste buitendiagonaalelement van A. 
Wegens symmetrie mogen we aannemen: p < q 


Ontwerp nu: 


T 
pq 5 pq: 


Deze matrix heeft natuurlijk dezelfde eigenwaarden als A. 


(H2,9) B = 0 


n n 
(#2,10) Opgave. Ga na dat z b?. = 5 Se (vgl. (42,5) ). 
daje 4 teder HJ 
9 kJ 


(42.11) Opgave. Ga na dat B weer symmetrisch is en dat 


D;;5 = Aij voor Ì * p, j * q 
: & . se 
„pd ° pj en: A 
= C - + 
„pp (PP a, ig KK 
ied s hand + „ 5 
vn win “pp daa B ê “pa 
aj * S 2pj * Cay voor j * p‚q 
bz b 
RS en ek 
qa * ® “pp * CS Apq * © äqa 
(42,12) Opgave. Ga na dat 
Se 28 (pp Ppa\ /° 8 bop _ Ppq 
a 8 
> 4 \ap Pay \ 5 Pap _ Paa 
(42,13) Opgave. Toon aan dat men c en s zo kan kiezen dat bra = bop z 0. 


(u2.14) Stel nu dat men c en s kiest als in (42,13). Dan geldt: 
tant vat stat 4 at 
bpp * Pqa * “pp pa “qa 
(vgl. (42,5) en (42.12). 


Omdat b;; a;j voor i f p,q geldt derhalve: 


2 e 2 2 
bij * as. + 2 a 


i=1 i 


Tan 3 


(42,15) De kwadraatsom der diagonaalelementen is aldus met 2 ana toegeno- 


men e 
De kwadraatsom der buitendiagonaalelementen moet dan met 2 ana 
afgenomen zijn (vgl. (42.10)), 


(42,16) Het prozes wordt nu herhaald met B i.p.v. A. 
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Doorgaans zullen nulgemaakte elementen bij de volgende stap weer 
* 0 worden. De kwadraatsom der buitendiagonaalelementen blijft 
echter zakken. 


(+7,17) ‘en hoopt natuurlijk dat de kwadraatsom der buitendiagonaalele- 
menten naar 0 gaat, waarna men eenvoudig van de diagonaal de ei- 
genwaarden kan aflezen. 

Dit zal zeker gebeuren als men telkens het absoluut grootste bui- 
tendiagonaalelement (vgl. (42.8) ) nul maakt. Zij immers 


N = n(n=-1) het aantal buitendiagonaalelementen, en zij S, de 
kwadraatsom der buitendiagonaalelementen na de k° transformatie- 
stap. 
Dan geldt max a: > S/N zodat 
AE k 
173 
Ss SS lied 0 
k+1 k ns 
p: 2,k+1 JN . 
Dus S44 * So (d-) en S4j ° 0 voor k +69, 


(42.18) Dit laatste zou overigens aanlviding kunnen zijn tot sombere ge= 
dachten over de convergentiesnelheid, maar dat valt mee, Men kan 
zelfs aantonen (moeilijk!) dat Jacobi superlineair convergeert. 


(42,19) In de praktijk blijken na circa 3n? transformatiestappen det bui- 
tendiagonaalelementen bij een rekenprecisie van 12 decimalen door- 
gaans niet meer significant van O te verschillen (empirisch). Aar:- 
gezien elke stap circa tn AV. kost is de totale hoeveelheid werk. 
dus omstreeks 12n?AV . Deze voorstelling van de hoeveelheid re- 
kenwerk is echter zeer bedriegelijk. Voor elke transformaticstanp 
moet men immers eerst het abs. grootste buitendiagonaalelement 
gaan bepalen. Bij het met de hand rekenen ziet men dat grootst. 
vlement op het oog, de rekenautomaat kan het grootste van Jn(n-i, 
elementen echter alleen deimn.v. Jn(n=-1)-1 aftrekkingen bepalen, 
dat is veel meer werk dan de bnAV. die vervolgens het uitvoeren 
van de transformatiestap kost! Vandaar dat men maar niet zoekt, 
maar gewoon achtereenvolypens Ajgs Ajzr eeen Aggeesesdgpsee sân, : 

nul maakt en dan weer van voren af aan begint (seriëel Jacobi 

proces). Mark op dat het nu helemaal niet meer duidelijk is dat 
het proces nog convergvert. 

Hen kan het echter wal bewijzen (niet eenvoudig). 
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(44.20) Corbato merkte op dat men om bij het oorspronkelijke Jacobi-pro- 


(ue. 


PG 


Naef 


ces het abs. grootste buitendiagonaalelement te bepalen toch nis: 
bij elke stap jn(n=1)=1 aftrekkingen hoeft uit te voeren. 

Per stap veranderen immers slechts 2 rijen en kolommen van de 
matrix. 

Door bij elke rij de plaats en de grootte van het abs, grootst. 
element van die rij te onthouden, heeft men per transformatiestaj. 
slechts circa 2n aftrekkingen nodig om het abs. grootste buiten- 
diagonaalelement te bepalen. Het is niet bekend of op deze wijze 
het echte Jacobiproees toch niet sneller convergeert dan het 


seriële proces. 


De numerieke stabiliteit van het Jacobiproces is goed. 
Stel dat men bij het rekenen in eindige precisie (dus na ongever 
3n? stappen voor het geval van 12 decimalen) uitkomt op. 
El 
4°8 8 \ 
En | 
\n G ei ) 
nz 
terwijl men bij exact rekenen had moeten krijgen: 
A N 





1 & 
RA 
IN 
© An 
Dan kan men aantonen dat geldt: 
\ 
n 
En Aj! À : 
\ iz1 sk < 2 = 
\ 5 108n E 
N EA. 
ì 


Voor ven 20 x 20 matrix betekent dit een verlies van U deeimale 
cijfers. bit is eehter cen absolute bovengrens, bij de bepaling 
waarvan men er rekening mee geiwouden heeft dat bij elke arithne 
tische operatie de maximale fout wordt gemaakt. Doordat de werk := 
lijke optredende fouten doorgaans (ondanks hun gedcetermincerdheid, 
ze zijn elke keer dat men hetzelfde proces draait, hetzelfde) cen 
stochastisch karakter lijken te hebben, treedt een soort statis= 
ike) effect op en wordt de werkelijke fout beter geschat door 

il n* &. 


Het Jacobiproces is niet het snelst bekende proces, wel het ven- 
voudigste, en daardoor voor niet te grote matrices zeer aantrek- 
kelijk. 


(42.23) 


(Hd. 24) 


(u2.25) 


(42,26) 


(u3.2) 
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Bovendien is het aantrekkelijk dat men meteen de eigenvectoren 
verkrijgt, nl. als de kolommen van het produkt der Oa (in de 
goede volgorde). 


Merk nog op dat men m.b.v. (36,12) kan zien hoe het # 0 zijn der 
buitendiagonaaluvlementen de evigenwaarden beïnvloedt en m.b.v. 


(37.28) hoe het mt de ei. envectoren staat. 


Opmerking. In principe kan men Jacobi toepassen op willekeurige 
nomaale matrices (vgl. (19.10) ), 

Een recent algorithme van Paardekooper transformeert stapsgewijs 
een willekeurige (niet noodzakclijk normale) matrix tot een ma- 
trix die bijna normaal is waarna men vervolgens Jacobi toepast 
om een (goede) benadering der eigenwaarden te vinden. 

M.H.C. Paardekooper = An Eigenvalue Algorithm based on Norm - 
Reducing Transformations. Diss Eindhoven 1969, 
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Methode van Houscholder - Givens. 


De methode van Householder = Givens behoort tot de tweede cate 
gorie van 638. 

Zij construcert in eindig vcel stappen een orthogonale matrix J 
die de reële symmetrische matrix A naar tridiagonaalvorm trans= 
formeert. 

De methode is oorspronkclijk ontwikkeld door Givens; later heef * 
Housecholder ven andere (snellere) manier aangegeven om de redue- 
tie tot tridiagonaalvorm tot stand te brengen. 


Bij Givens wordt O stapsgewijs opgebouwd uit draaiingen zoals 
in (k2,6). Icdere draaiing is bedoeld om 2 symmetrisch gelegen 
matrixclementen buiten de tridiagonale band nul te maken, waar- 
bij elvmenten die eenmaal nulgemaakt zijn nul blijven gedurende 
het verdere verloop van het proces. 


(u3,3) 


(U3.4) 


(H3,6) 


(k 3.7) 


De methode gaat als volgt: 


=ecerst worden as en a13 nulgemaakt door een rotatie in het 


‚…) 
| 
„dan worden ba, en by4 nulgemaakt door een rotatie in het 


(eas e3)-vlak,. Resulterende matrix B = (b, 


(ea , va )d=vlak. 


„cete, 


Opgave. Ga na dat nulgemaakte clementen inderdaad nul blijven. 


Bij Housecholder wordt O stapsgewijs opgebouwd uit Houscholder- 
transformaties Hz I= 2w.wl waarin w een vector met 2enorm 1 
is (vele (43,3) Js 


In de eerste stap Houscholder bipaalt men Hs (dus cigenlijk w) 
zó dat Azie Age vers Ang PR Age Agur wees Aan alle tegelijk 
nul worden. Resulterende matrix: Bz (bij): 
In de tweed, stap Houscholder bepaalt men HK, zo dat de 2lementer 


bijz» eg bz 
bjjs wees Dog Sn Djgs vee Dj, NUL blijven). Etc. 


en bou» bed bon nul worden (en de elementen 


De benodigde rekentijd voor de reductie tot tridiagonaalvorm is 


beperkt : n° AV bij Givens en an AV bij Houscholder. 


Beide algorithmen zijn zeer stabiel. 
Noteren we met Aj resp. A de (naar grootte gerangschikte) eigen- 
waarden van de oorspronkelijke en van de in eindige precisie berekende 
tridiagonaalmatrix dan, geldt 

Eat? |} 
1 Ì i 3/ 


Ä zin 2 E bij Givens 





E42 
k 8 


n | S on? € bij Houscholder 
(vergelijk ook (42,21) ). 


Zij T de verkregen tridiagonaalmatrix en stel 
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8 tet een A 
A4 b4 8} 
b4 ä, \ 
= is, I 4 
ne | | 
T = | i k: ik is Ee 0 
\ 1 f En End nd 
EE: Ne en Dr-1/ 
0-0 boat än 


De waarde van Y(T;A) = det (T- AI) is nu heel goedkoop te bee 
palen via cen recurrente betrekking. 


(43.8) Definiter nl. 


| 
Aj d. 
b iN | 

jn Ne ss 

p;(Â) = 1 72 % | (i 21) 
N 
biea | 
bg Ak | 
Met P_, 0) z= 0 » Po) = 1 krijgt men: 
( = a … he 3 7 
(43,9) P;(A) = (a, = A) p;_g OÙ = big Pig) Cis 1} 


43.1c) Opgave. Verifiter (53,8) en ga na dat VT$A) z PA). 


(63,11) Evenwel, in plaats van nu bijv. Koorden = Newton toe te passen 
om de nulpunten van Y(T;A) te bepalen, gebruikt men de rij 
{p; A) } op cen andere manier. 
Onder aanname dat alle bd, %* Olwat betekent b; z 02) vormt de rij 
Pps Pj°***+s Pp een zgn. Sturmerij van polynomen, d.w.z. geen van 
de p; is = 0, en in ven nulpunt @ van p; (i > 1) geldt 
sEn(p;_4(9) ) = =sgn(p;,,(9) ) #0, Zie verder gaa 


W3,12) Opgave. Verifiter dit laatste m.b.v. (43,9) 


43.13) Met behulp van de Sturm-rij kan men ven heel bruikbaar procédé 
voor de bepaling van de evigenwaarden van T ontwerpen. 
Verdere bijzonderheden vindt men bv. bij Wilkinson, Ortega en 


Givens. 


(H3.18) Het beprlen van de cigenveetor:n van ven tridiagonaalmatrix is 
cen vervelende zaak. Hiervoor is nl, geen stabiele algorithme be- 
kend. (vgl. Wilkinson p. 315 c.v.). Soms gebruikt men Wielandt 
iteratic,. 





(uut) Ala alls Ho dn UAC, me Aan es cl a saf, 


(44,2) 


(44,3) 


(4u 4) 
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TENTAMEN LINEAIRE ÄLGEBRA 
(Numerieke Analyse II) 


ZATERDAG 8 MAART 1975 3.30, 30 


es LEES ONDERSTAANDE INSTRUKTIES GOED DOOR 

ie ZET OP ELK INGELEVERD VEL UW NAAM 

hi ZET OP HET DUBBELE FOLIOVEL ALLEEN UW NAAM+ADRES 
DIT VEL MAG VERDER NIET BESCHREVEN WORDEN EN 
DIENT OM UW OPGAVEN IN TE LEVEREN 

a U MAG VOORGAANDE ONDERDELEN GEBRUIKEN, OOK ALS 
ZE NIET BEWEZEN ZIJN, 


en 


[no 


We willen bewijzen dat voor alle reële nxn matrices 
A geldt: 
WA, = Al > rang(A) <1. 
Zoals bekend bestaan er orthogonale matrices U en V en 
T 


een niet negatieve diagonaalmatrix XZ zodat A=UDV". 


Ga nu achtereenvolgens na: 


(a) Zl, = All2, LZ = Al 
(b) rang(X) = rang(A) (bedenk dat rang(A) = dim(AR* )) 


Ce) All, = LAI > pang(A) S 1. 


dû 4 4 
Zij A = 1 dU A 
OE 


(a) Geef afschattingen voor de eigenwaarden van A mbv. 
de stellingen van Gershgorin. 

(b) Bepaal een diagonaal matrix D zodanig dat de stellingen 
van Gershgorin toegepast op DAD” Î aantonen dat er een 
eigenwaarde tussen 9.7 en 10.3 ligt. 

(ec) Laat zien dat O0 een eigenwaarde is van A, en toon 
vervolgens aan, zonder de karakteristieke veelterm op te 
stellen, dat de derde eigenwaarde tussen 1.7 en 2,3 ligt. 

(d) Men kan A ook beschouwen als een symmetrisch gepertur= 
beerde van diag (10, 1, 1). Dit geeft de mogelijkheid de 
eigenwaarden te schatten met behulp van de 2-norm van de 
stoormatrix (die hiertoe exact bepaald moet worden). 


Doe dit. 


joo 
e 


Zij A een mxm matrix, b een gegeven vector. 

We wensen de oplossing van Ax = b benaderd te bepalen. 
Als norm kiezen we de supnorm,. 

(a) Zij Xx een benadering van x., Bewijs nu dat 

Ixil SHA UPG met P(X) = Ak-b. Ga ook na dat indien 
LAS eG SIX er geldt: 


x= HAT ÀL 50 
En , . 
WX AT Ie GO 


(b) Stel nu dat xo de oplossing van het stelsel Bx, = b is 
(B non-singulier!), waarbij IL 1-8” tal < 1 is. 
Toon aan dat de rij vectoren {x_} met 


Ge), = (I-B ÎA)x, +5 Ib (n>1) naar de exacte oplossing van 


„1 


Ax = b convergeert (afgezien van afrond fouten). 


(e) Zij B de matrix E Demel a: 
ee a ® 
ER 
hi 
Ò Pel 
‚k+T <1 
a 1 
m 


(d) Geef een expliciete uitdrukking voor B°Î en toon aan dat 


(Bp = 1+max heil 
1Sj Sn 
(e) Zij nu A = B+E met IENS e en e(1+max |a;|) < 0.1 
1SjSn } et 
Geef een zo scherp mogelijke afschatting voor ||I-B "All en 


VAT Ä, en laat zien dat in elk geval Ir-BÎAl < 0.1 


(£) Stel B'‘b is exact te bepalen; noteer LB Îolt = ce. 


9 10 
an 
Toon aan Toe S xl < Te 


(g) 


Ch) 


Ci) 


(5) 


Zij P een willekeurige computer representeerbare mXxm matrix 
en Xx een dito vector, Laat zien dat de norm van de fout bij 
het berekenen van Px hoogstens miiPlillxllE is als & de max. 
relatieve machine afrondfout is en we hogere machten van E 


verwaarlozen. 


1 1 


Neem nu aan dat de matrix I-B “A en de vector B “b exact 


bekend zijn en door de computer representeerbaar, en laat de 
computer de recursie (#) ûreuoeen met deze exacte ie 
Opgeleverd wordt dan een rij lek, 

Toon aan 


Wx_-[CT-B'ÎA)x_4+B DI < [o,1lmti)llx,_slltelE, 


-1 
als we hogere machten van £ verwaarlozen, 


Zij Ö = KE Toon dan mbv, voorgaande resultaten en (+) aan 


1 5 
Hell S 0,18 UL + 5 (m+10)e & 


„1 
en laat hiermee zien dat Wôl S El (m+10)e E 


Toon nu aan dat voor n groot genoeg zeker geldt 


EN 


0 


Wx =xl 5 (m+10)e E. 


= 8 


Hal 
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